1 SETUP

1.1 Open libraries

library(gtsummary)
library(dplyr)
library(tidyverse)
library(survival)
library(survminer)
library(lubridate)
library(readxl)
library(janitor)
library(ggplot2)
library(tidycmprsk)
library(ggsurvfit)
library(tidymodels) 
library(sjPlot)
library(patchwork)
library(bstfun)
library(ComplexUpset)

library(writexl)
library(rms)

library(timeROC)
library(rms)
library(pec)
library(prodlim)
library(riskRegression)

library(plotly)
library(ggtext)

library(WeightIt)
library(cobalt)
library(survey)

2 DATA WRANGLING

2.1 Read in MM datasets

db_Kansas = read_excel(path = "/Users/aportugu/Desktop/Elranatamab/Analysis/Data/bsAb Consortium Data Template 2.11.2025 - due 5-1-25 (de-identified).xlsx", sheet = "bsAb_DATA")%>%
  filter(`bsAb name` %in% c("1","3"))%>%
  clean_names()%>%
  filter(!is.na(age_at_bs_ab_initiation)) %>%
  mutate(
    patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate = as.character(patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate),
    baseline_anc = as.double(baseline_anc),
    baseline_alc = as.double(baseline_alc),
    alc_day_8 = as.double(alc_day_8),
    alc_day_15 = as.double(alc_day_15),
    alc_day_22 = as.double(alc_day_22),
    alc_c2d1 = as.double(alc_c2d1),
    anc_nadir_value_during_first_6_cycles = as.double(anc_nadir_value_during_first_6_cycles),
    cr_cl_prior_to_step_up_dose_number_1 = as.double(cr_cl_prior_to_step_up_dose_number_1)
  )

# Write to Excel
#write_xlsx(db_Kansas, "Kansas_elra_data.xlsx")

db_CC = read_excel(path = "/Users/aportugu/Desktop/Elranatamab/Analysis/Data/CC bsAb Consortium Submission CCF 6-11-25.xlsx", sheet = "bsAb_DATA")%>%
  filter(`bsAb name` %in% c("1","3"))%>%
  clean_names()%>%
  mutate(
    patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate = as.character(patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate),
    baseline_anc = as.double(baseline_anc),
    baseline_alc = as.double(baseline_alc),
    alc_day_8 = as.double(alc_day_8),
    alc_day_15 = as.double(alc_day_15),
    alc_day_22 = as.double(alc_day_22),
    alc_c2d1 = as.double(alc_c2d1),
    anc_nadir_value_during_first_6_cycles = as.double(anc_nadir_value_during_first_6_cycles),
    cr_cl_prior_to_step_up_dose_number_1 = as.double(cr_cl_prior_to_step_up_dose_number_1)
  )

# Write to Excel
#write_xlsx(db_CC, "CC_elra_data.xlsx")

db_HUMC = read_excel(path = "/Users/aportugu/Desktop/Elranatamab/Analysis/Data/BsAb Consortium Data HUMC.xlsx", sheet = "bsAb_DATA")%>%
  filter(`bsAb name` %in% c("1","3"))%>%
  clean_names()%>%
  mutate(
    patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate = as.character(patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate),
    baseline_anc = as.double(baseline_anc),
    baseline_alc = as.double(baseline_alc),
    alc_day_8 = as.double(alc_day_8),
    alc_day_15 = as.double(alc_day_15),
    alc_day_22 = as.double(alc_day_22),
    alc_c2d1 = as.double(alc_c2d1),
    anc_nadir_value_during_first_6_cycles = as.double(anc_nadir_value_during_first_6_cycles),
    cr_cl_prior_to_step_up_dose_number_1 = as.double(cr_cl_prior_to_step_up_dose_number_1)
  )

# Write to Excel
#write_xlsx(db_HUMC, "HUMC_elra_data.xlsx")

db_MCC = read_excel(path = "/Users/aportugu/Desktop/Elranatamab/Analysis/Data/bsAb Consortium Data MCC-deidentified.xlsx", sheet = "bsAb_DATA")%>%
  filter(`bsAb name` %in% c("1","3"))%>%
  clean_names()%>%
  mutate(
    patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate = as.character(patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate),
    anc_nadir_value_during_first_6_cycles = as.double(anc_nadir_value_during_first_6_cycles),
    baseline_anc = as.double(baseline_anc),
    baseline_alc = as.double(baseline_alc),
    cr_cl_prior_to_step_up_dose_number_1 = as.double(cr_cl_prior_to_step_up_dose_number_1)
  )


db_Huntsman = read_excel(path = "/Users/aportugu/Desktop/Elranatamab/Analysis/Data/Utah_Huntsman_bsAb Consortium Data FINAL_5.4.25.xlsx", sheet = "bsAb_DATA")%>%
  filter(`bsAb name` %in% c("1","3"))%>%
  clean_names()%>%
  mutate(
    patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate = as.character(patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate),
    baseline_anc = as.double(baseline_anc)/1000,
    baseline_alc = as.double(baseline_alc)/1000,
    anc_nadir_value_during_first_6_cycles = as.double(anc_nadir_value_during_first_6_cycles),
    cr_cl_prior_to_step_up_dose_number_1 = as.double(cr_cl_prior_to_step_up_dose_number_1)
  )

# Write to Excel
#write_xlsx(db_Huntsman, "Huntsman_elra_data.xlsx")


db_FHCC = read_excel(path = "/Users/aportugu/Desktop/Elranatamab/Analysis/Data/FHCC Bispecifics (5.5.2025).xlsx", sheet = "bsAb_DATA")%>%
  filter(`bsAb name` %in% c("1","3"))%>%
  clean_names()%>%
  mutate(
    patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate = as.character(patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate),
    date_of_infection_80 = ymd(date_of_infection_80),
    baseline_anc = as.double(baseline_anc),
    cr_cl_prior_to_step_up_dose_number_1 = as.double(cr_cl_prior_to_step_up_dose_number_1)
  )


db_UTSW = read_excel(path = "/Users/aportugu/Desktop/Elranatamab/Analysis/Data/bsAb Consortium Data Template 2.11.2025 1.xlsx", sheet = "bsAb_DATA")%>%
  filter(`bsAb name` %in% c("1","3")) %>%
  clean_names()%>%
  mutate(
    patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate = as.character(patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate),
    baseline_anc = as.double(baseline_anc),
    baseline_alc = as.double(baseline_alc),
    anc_nadir_value_during_first_6_cycles = as.double(anc_nadir_value_during_first_6_cycles)
  )

# Write to Excel
#write_xlsx(db_UTSW, "UTSW_elra_data.xlsx")


db_MDACC = read_excel(path = "/Users/aportugu/Desktop/Elranatamab/Analysis/Data/MDACC BITE ASH 2025.xlsx", sheet = "bsAb_DATA")%>%
  filter(`bsAb name` %in% c("1","3"))%>%
  clean_names()%>%
  mutate(
    patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate = as.character(patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate),
    baseline_anc = as.double(baseline_anc),
    baseline_alc = as.double(baseline_alc),
    anc_nadir_value_during_first_6_cycles = as.double(anc_nadir_value_during_first_6_cycles),
    cr_cl_prior_to_step_up_dose_number_1 = as.double(cr_cl_prior_to_step_up_dose_number_1),
    total_number_paraskelatal_plasmacytomas = as.double(total_number_paraskelatal_plasmacytomas)
  )

# Write to Excel
#write_xlsx(db_MDACC, "MDACC_elra_data.xlsx")


db_Roswell = read_excel(path = "/Users/aportugu/Desktop/Elranatamab/Analysis/Data/bsAb Consortium Data Template 5.8.2025_NN.xlsx", sheet = "bsAb_DATA")%>%
  filter(`bsAb name` %in% c("1","3"))%>%
  clean_names()%>%
  mutate(
    patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate = as.character(patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate),
    baseline_anc = as.double(baseline_anc),
    baseline_alc = as.double(baseline_alc),
    anc_nadir_value_during_first_6_cycles = as.double(anc_nadir_value_during_first_6_cycles),
    cr_cl_prior_to_step_up_dose_number_1 = as.double(cr_cl_prior_to_step_up_dose_number_1),
    total_number_paraskelatal_plasmacytomas = as.double(total_number_paraskelatal_plasmacytomas)
  )

# Write to Excel
#write_xlsx(db_Roswell, "Roswell_elra_data.xlsx")

# Combine all into a single dataframe
combined_db <- bind_rows(
  list(
    HUMC = db_HUMC,
    MCC = db_MCC,
    Huntsman = db_Huntsman,
    FHCC = db_FHCC,
    UTSW = db_UTSW,
    MDACC = db_MDACC,
    Roswell = db_Roswell,
    CC = db_CC,
    Kansas = db_Kansas
  ),
  .id = "site"
) %>%
  mutate(
    prior_bcma_directed_therapy_yes_no = ifelse(!is.na(most_recent_bcma_agent_name),1,0),
    bs_ab_name = case_when(
      bs_ab_name == 1 ~ "tec",
      bs_ab_name == 3 ~ "elra",
      .default = NULL
    ),
    max_crs_grade_0_5 = ifelse(is.na(max_crs_grade_0_5), 0,max_crs_grade_0_5), ## Set NAs to 0
    max_icans_grade = ifelse(is.na(max_icans_grade), 0, max_icans_grade ) ## Set NAs to 0
  )

# Convert all datetime columns to Date
combined_db_clean <- combined_db %>%
  mutate(across(where(lubridate::is.POSIXt), as.Date))

# Write to Excel
write_xlsx(combined_db_clean, "elra_tec_data.xlsx")


# Write to Excel
# write_xlsx(combined_db_clean %>% filter(prior_bcma_directed_therapy_yes_no != prior_bcma_directed_therapy), "elra_data_discordant.xlsx")

## Several institutions do not have any cases of elranatamab

# db_UABMC = read_excel(path = "Data/bsAb Consortium Data 5.8.xlsx", sheet = "bsAb_DATA")%>%
#   filter(`bsAb name` == "3")%>%
#   clean_names()

# db_Stanford = read_excel(path = "Data/De-identified Stanford bsAb Consortium Data Template 4.30.2025.xlsx", sheet = "bsAb_DATA")%>%
#   filter(`bsAb name` == "3")%>%
#   clean_names()

# db_LCI = read_excel(path = "Data/bsAb Consortium Data LCI.xlsx", sheet = "bsAb_DATA")%>%
#   filter(`bsAb name` == "3")%>%
#   clean_names()

# db_COH = read_excel(path = "Data/bsAb Consortium Data Template 2.11.2025_COHv2.xlsx", sheet = "bsAb_DATA")%>%
#   filter(`bsAb name` == "3")%>%
#   clean_names()

# db_MUSC = read_excel(path = "Data/USMM bsAb Consortium Data Template 2.11.2025.xlsx", sheet = "bsAb_DATA")%>%
#   filter(`bsAb name` == "3")%>%
#   clean_names()

# db_Mayo = read_excel(path = "Data/Bispecific_consortium_20250521 NoPHI.xlsx", sheet = "bsAb_DATA") %>%
#   filter(`bsAb name` == 3)%>%
#   clean_names()

2.2 Calculate CAR-HEMATOTOX score and risk level

combined_db <- combined_db %>% 
  mutate(
    # baseline_anc = as.numeric(baseline_anc),
    # baseline_hgb,
    # baseline_pl_ts,
    # baseline_crp_prior_to_bs_ab_initiation,
    # baseline_ferritin_prior_to_bs_ab_initiation,
    HT_Score = 
      (baseline_anc <= 1.2) +
      (baseline_hgb <= 9.0) +
      (baseline_pl_ts >= 76 & baseline_pl_ts <= 175) +
      (baseline_crp_prior_to_bs_ab_initiation >= 3.0) +
      (baseline_ferritin_prior_to_bs_ab_initiation >= 650 & baseline_ferritin_prior_to_bs_ab_initiation <= 2000) +
      (baseline_pl_ts <= 75) + ## Add an extra point for Plt <=75
      (baseline_ferritin_prior_to_bs_ab_initiation > 2000), ## Add an extra point for ferritin >2000
    
    # Classify the risk based on the score
    HT_Risk = factor(
      case_when(
        HT_Score >= 2 ~ "HThigh",
        baseline_pl_ts <=75 ~ "HThigh", ## Equivalent to at least 2 points
        baseline_ferritin_prior_to_bs_ab_initiation > 2000 ~ "HThigh", ## Equivalent to at least 2 points
        (baseline_anc <= 1.2) + (baseline_hgb <= 9.0) + (baseline_pl_ts >= 76 & baseline_pl_ts <= 175) >=2 ~ "HThigh",
        HT_Score <2 ~ "HTlow",
        .default = "NA"
      ),
      levels = c("HTlow", "HThigh")
    )
)
write_xlsx(combined_db %>% select(site, patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate, baseline_anc, baseline_hgb, baseline_pl_ts, baseline_crp_prior_to_bs_ab_initiation, baseline_ferritin_prior_to_bs_ab_initiation, HT_Score, HT_Risk), "HT_score.xlsx")

2.3 Identify severe infections

combined_db <- combined_db %>% 
  mutate(
    severe_infection = ifelse(infection_severity_73 %in% 1 | infection_severity_76 %in% 1 | infection_severity_79 %in% 1| infection_severity_82 %in% 1,1,0),
    first_severe_infection_date = case_when(
      infection_severity_73 %in% 1 ~ date_of_infection_74,
      infection_severity_76 %in% 1 ~ date_of_infection_77,
      infection_severity_79 %in% 1 ~ date_of_infection_80,
      infection_severity_82 %in% 1 ~ date_of_infection_83,
      .default = NULL
    ),
    first_severe_infection_type = case_when(
      infection_severity_73 %in% 1 ~ number_infection_1_bacterial_viral_fungal,
      infection_severity_76 %in% 1 ~ number_infection_2_bacterial_viral_fungal,
      infection_severity_79 %in% 1 ~ number_infection_3_bacterial_viral_fungal,
      infection_severity_82 %in% 1 ~ number_infection_4_bacterial_viral_fungal,
      .default = NULL
    )
    
    
  )

3 PROPENSITY SCORES

3.1 Regression analysis for treatment allocation

theme_gtsummary_compact()

variables <- c("elra", "Age", "Female sex", "ECOG 2+", "EMD", "Hgb", "LDH", "Prior BCMA")

data1 <- combined_db %>% 
  drop_na(age_at_bs_ab_initiation, ecog_ps_at_bs_ab_initiation_0_5, just_prior_to_bs_ab_initiation_is_there_true_emd_if_none_move_to_column_gk, baseline_hgb, baseline_ldh_prior_to_bs_ab_initiation) %>%
  transmute(
    Age = age_at_bs_ab_initiation,
    "Female sex" = sex_0_male_1_female,
    elra = ifelse(bs_ab_name == "elra", 1, 0),
    "ECOG 2+" = ifelse(ecog_ps_at_bs_ab_initiation_0_5 >=2,1,0),
    EMD = ifelse(just_prior_to_bs_ab_initiation_is_there_true_emd_if_none_move_to_column_gk == 2,1,0),
    Hgb = as.double(baseline_hgb),
    LDH = as.double(baseline_ldh_prior_to_bs_ab_initiation),
    "Prior BCMA" = prior_bcma_directed_therapy_yes_no
  )

uv_tab <- tbl_uvregression(
    data1[c(variables)],
    method = glm,
    y = "elra",
    method.args = list(family = binomial),
    exponentiate = TRUE
  ) %>%
  sort_p()


mv_tab<-glm(elra ~ Age + `Female sex` + `ECOG 2+` + EMD + Hgb + LDH + `Prior BCMA`, data=data1, family=binomial) %>%
  tbl_regression(exponentiate=TRUE)

tbl_merge(list(uv_tab, mv_tab), tab_spanner = c("**Univariate**", "**Multivariable**"))
Characteristic
Univariate
Multivariable
N OR1 95% CI1 p-value OR1 95% CI1 p-value
Hgb 371 1.12 1.01, 1.25 0.029 1.14 1.02, 1.28 0.026
EMD 371 1.69 0.95, 2.96 0.069 1.71 0.96, 3.04 0.068
Age 371 1.02 1.00, 1.04 0.087 1.01 0.99, 1.04 0.3
Female sex 371 1.43 0.93, 2.21 0.11 1.43 0.92, 2.24 0.11
Prior BCMA 371 0.76 0.49, 1.17 0.2 0.78 0.50, 1.23 0.3
ECOG 2+ 371 1.19 0.75, 1.88 0.4 1.27 0.78, 2.07 0.3
LDH 371 1.00 1.00, 1.00 0.8 1.00 1.00, 1.00 >0.9
1 OR = Odds Ratio, CI = Confidence Interval

3.2 IPTW

IPTW.data <- combined_db %>% 
  drop_na(age_at_bs_ab_initiation, ecog_ps_at_bs_ab_initiation_0_5, just_prior_to_bs_ab_initiation_is_there_true_emd_if_none_move_to_column_gk, baseline_hgb, baseline_ldh_prior_to_bs_ab_initiation) %>%
  mutate(
    Age = age_at_bs_ab_initiation,
    "Female sex" = sex_0_male_1_female,
    bsAb = bs_ab_name,
    "Received full dose" = ifelse(!is.na(date_of_first_bs_ab_full_dose),1,0),
    "ECOG" = case_when(
      ecog_ps_at_bs_ab_initiation_0_5 %in% c("0","1") ~"0-1",
      ecog_ps_at_bs_ab_initiation_0_5 ==2 ~"2",
      ecog_ps_at_bs_ab_initiation_0_5 >=3 ~"3+",
      .default = NULL),
    EMD = ifelse(just_prior_to_bs_ab_initiation_is_there_true_emd_if_none_move_to_column_gk == 2,1,0),
    Hgb = as.double(baseline_hgb),
    LDH = as.double(baseline_ldh_prior_to_bs_ab_initiation),
    "Prior BCMA" = prior_bcma_directed_therapy_yes_no
  )

## Use the WeightIt function to generate propensity scores and weights
W.out <- weightit(bsAb ~ Age + `ECOG` + EMD + Hgb + LDH + `Prior BCMA`,
                  data=IPTW.data, 
                  method = "glm",    ## Methods: glm, gbm, cbps, npcbps, ebal, optweight, super, bart, energy
                  estimand = "ATE",
                  stabilize = TRUE)

bal.tab(W.out, stats = c("m", "v"), thresholds = c(m=0.25), binary = "std")
## Balance Measures
##                Type Diff.Adj     M.Threshold V.Ratio.Adj
## prop.score Distance   0.0050 Balanced, <0.25      0.9733
## Age         Contin.   0.0090 Balanced, <0.25      0.9595
## ECOG_0-1     Binary  -0.0121 Balanced, <0.25           .
## ECOG_2       Binary   0.0152 Balanced, <0.25           .
## ECOG_3+      Binary  -0.0027 Balanced, <0.25           .
## EMD          Binary  -0.0034 Balanced, <0.25           .
## Hgb         Contin.  -0.0186 Balanced, <0.25      1.0974
## LDH         Contin.  -0.0276 Balanced, <0.25      1.0961
## Prior BCMA   Binary  -0.0119 Balanced, <0.25           .
## 
## Balance tally for mean differences
##                     count
## Balanced, <0.25         9
## Not Balanced, >0.25     0
## 
## Variable with the greatest mean difference
##  Variable Diff.Adj     M.Threshold
##       LDH  -0.0276 Balanced, <0.25
## 
## Effective sample sizes
##              elra    tec
## Unadjusted 123.   248.  
## Adjusted   113.38 242.39
summary(W.out)
##                   Summary of weights
## 
## - Weight ranges:
## 
##         Min                                  Max
## elra 0.4745 |---------------------------| 2.3134
## tec  0.7811      |-------------|          1.7418
## 
## - Units with the 5 most extreme weights by group:
##                                         
##          123    102     81     91      9
##  elra 1.5458 1.5896  1.653 1.7607 2.3134
##          169     71    274    258    271
##   tec 1.4257 1.4453 1.5215 1.6783 1.7418
## 
## - Weight statistics:
## 
##      Coef of Var   MAD Entropy # Zeros
## elra       0.292 0.222   0.040       0
## tec        0.152 0.115   0.011       0
## 
## - Effective Sample Sizes:
## 
##              elra    tec
## Unweighted 123.   248.  
## Weighted   113.38 242.39
love_plot <- love.plot(W.out,
          binary = "std",
          line = TRUE, 
          abs = TRUE,
          thresholds = c(m=0.2),
          sample.names = c("Unweighted", "IPTW"),
          title = "A)",
          var.order = "unadjusted",
            var.names = c(
              prop.score = "Distance"
          ))+ theme(plot.title = element_text(hjust = 0))
love_plot

weight_df = data.frame(W.out$weights, W.out$treat)

IPTW.data$weights <- W.out$weights
  
IPTW.data$ps <- W.out$ps


write.csv(data1, "IPTW.dataset.csv")

3.3 Scatterplot of weights

Weights_plot <- ggplot(
  IPTW.data, aes(x=ps, y = weights, color = bsAb)) + 
  geom_point(shape = 16) +
  theme_classic() + 
  scale_y_continuous(breaks = seq(0,4 , by = 0.5)) +
  labs(x = "Propensity score", y = "Weights") + 
  ggtitle("B)")
Weights_plot

4 TABLES

4.1 Site contributions

theme_gtsummary_compact()

combined_db %>%
  transmute(
    "Site" = site,
    bs_ab_name,
  ) %>%
  tbl_summary(
    by = bs_ab_name,
    missing = "no",
    sort = list(all_categorical() ~ "frequency"),
    statistic = list(
      all_continuous() ~ "{median} ({p25}, {p75})",
      all_categorical() ~ "{n} ({p}%)"
    )
  ) %>%
  bold_labels()
Characteristic elra
N = 130
1
tec
N = 277
1
Site

    MCC 53 (41%) 72 (26%)
    MDACC 12 (9.2%) 57 (21%)
    CC 20 (15%) 47 (17%)
    UTSW 10 (7.7%) 40 (14%)
    FHCC 7 (5.4%) 24 (8.7%)
    Huntsman 9 (6.9%) 22 (7.9%)
    Kansas 4 (3.1%) 13 (4.7%)
    HUMC 10 (7.7%) 1 (0.4%)
    Roswell 5 (3.8%) 1 (0.4%)
1 n (%)

4.2 Unweighted and IPTW baseline characteristics

theme_gtsummary_compact()

data1 = IPTW.data %>%
  transmute(
    "Age (years)" = age_at_bs_ab_initiation,
    "Female sex" = sex_0_male_1_female == 1,
    "Race" = factor(case_when(
      is.na(race_white_0_black_1_aapi_2_american_indian_3_other_4_unk_5) ~ NA,
      race_white_0_black_1_aapi_2_american_indian_3_other_4_unk_5 == 0 ~ "White",
      race_white_0_black_1_aapi_2_american_indian_3_other_4_unk_5 == 1 ~ "Black",
      race_white_0_black_1_aapi_2_american_indian_3_other_4_unk_5 == 2 ~ "AAPI",
      race_white_0_black_1_aapi_2_american_indian_3_other_4_unk_5 == 3 ~ "American Indian",
      race_white_0_black_1_aapi_2_american_indian_3_other_4_unk_5 == 4 ~ "Other",
      race_white_0_black_1_aapi_2_american_indian_3_other_4_unk_5 == 5 ~ NA
    ), levels = c("White","Black","AAPI","American Indian","Other",NA)),
    #ecog_ps_at_bs_ab_initiation_0_5,
    "ECOG" = case_when(
      ecog_ps_at_bs_ab_initiation_0_5 == 0 ~ "0",
      ecog_ps_at_bs_ab_initiation_0_5 == 1 ~ "1",
      ecog_ps_at_bs_ab_initiation_0_5 == 2 ~ "2",
      ecog_ps_at_bs_ab_initiation_0_5 >= 3 ~ "3+"
    ),
    "LDH (U/L)" = as.double(baseline_ldh_prior_to_bs_ab_initiation),
    "Hgb (g/dL)" = baseline_hgb,
    "Baseline ferritin (ng/mL)" = as.double(baseline_ferritin_prior_to_bs_ab_initiation),
    "Heavy chain" = factor(case_when(
      myeloma_involved_heavy_chain_ig_g_ig_a_ig_m_ig_d_ig_e_or_blank == "IgG" ~ "IgG",
      myeloma_involved_heavy_chain_ig_g_ig_a_ig_m_ig_d_ig_e_or_blank == "IgM" ~ "IgM",
      myeloma_involved_heavy_chain_ig_g_ig_a_ig_m_ig_d_ig_e_or_blank == "IgA" ~ "IgA",
      myeloma_involved_heavy_chain_ig_g_ig_a_ig_m_ig_d_ig_e_or_blank == "IgD" ~ "IgD",
      .default = "None"
      ), levels = c("IgG", "IgA","IgM","IgD","None")),
    "True EMD" = just_prior_to_bs_ab_initiation_is_there_true_emd_if_none_move_to_column_gk == 2,
    "≥1 HRCA" = ifelse(deletion_17p_prior_to_bs_ab_0_no_1_yes==1 | t_4_14_prior_to_bs_ab_0_no_1_yes == 1 | t_14_16_prior_to_bs_ab_0_no_1_yes == 1 | amp_1_q_only_0_no_1_yes == 1 | gain_amp1q21_prior_to_bs_ab_0_no_1_yes == 1,1,0),
    "Prior LOTs" = number_of_prior_lines_of_therapy_before_bs_ab,
    "Prior ASCT" = prior_auto_sct_no_0_yes_1,
    "Prior CAR-T" = prior_cart_0_no_1_yes,
    "Triple refractory" = triple_class_refractory_i_mi_d_pi_anti_cd38 == 1,
    "Penta refractory" = penta_refractory_2_p_is_2_i_mi_ds_anti_cd38 == 1,
    "Prior BCMA-directed therapy" = !is.na(most_recent_bcma_agent_name),
    "Trial eligible" = eligible_for_corresponding_bs_ab_trial_on_first_dose_use_trial_elig_sheet_to_help_you == 1,
    "Received full dose" = ifelse(!is.na(date_of_first_bs_ab_full_dose),1,0),
    bs_ab_name,
    weights
  )

table_unweighted <- data1 %>% select(-weights) %>%
  tbl_summary(
    by = bs_ab_name,
    missing = "no",
    #sort = list(all_categorical() ~ "frequency"),
    statistic = list(
      all_continuous() ~ "{median} ({p25}, {p75})",
      all_categorical() ~ "{n} ({p}%)",
      "Age (years)" ~ "{median} ({min} to {max})"
    )
  ) %>%
  add_p() %>%
  bold_labels()

table_IPTW <- survey::svydesign(
    id = ~0,
    weights = ~weights,
    data = data1,
    fpc = NULL
  ) %>%
  tbl_svysummary(
    by = bs_ab_name,
    missing = "no",
    statistic = list(
      all_continuous() ~ "{median} ({p25}, {p75})",
      all_categorical() ~ "{n} ({p}%)"),
    include = c("Age (years)", "Female sex","Race","ECOG","LDH (U/L)","Hgb (g/dL)","Baseline ferritin (ng/mL)", "Heavy chain","True EMD","≥1 HRCA","Prior LOTs","Prior ASCT","Prior CAR-T","Triple refractory","Penta refractory","Trial eligible","Prior BCMA-directed therapy", "Received full dose")
  ) %>%
  add_p() %>%
    bold_labels()
  
tbl_merge(list(table_IPTW, table_unweighted), tab_spanner = c("**IPTW**", "**Unweighted**"))
Characteristic
IPTW
Unweighted
elra
N = 123
1
tec
N = 248
1
p-value2 elra
N = 123
3
tec
N = 248
3
p-value4
Age (years) 70 (64, 76) 69 (62, 76) 0.9 71 (39 to 95) 68 (35 to 88) 0.060
Female sex 66 (54%) 114 (46%) 0.2 68 (55%) 115 (46%) 0.11
Race

0.2

0.15
    White 94 (79%) 187 (76%)
95 (81%) 187 (76%)
    Black 17 (14%) 50 (20%)
16 (14%) 51 (21%)
    AAPI 3 (2.5%) 6 (2.3%)
3 (2.5%) 5 (2.0%)
    American Indian 2 (1.8%) 0 (0%)
2 (1.7%) 0 (0%)
    Other 2 (2.0%) 4 (1.8%)
2 (1.7%) 4 (1.6%)
ECOG

>0.9

0.5
    0 23 (18%) 49 (20%)
21 (17%) 50 (20%)
    1 62 (50%) 119 (48%)
59 (48%) 121 (49%)
    2 28 (23%) 59 (24%)
29 (24%) 60 (24%)
    3+ 10 (8.2%) 20 (8.2%)
14 (11%) 17 (6.9%)
LDH (U/L) 223 (178, 280) 208 (167, 290) 0.4 222 (175, 277) 208 (167, 291) 0.5
Hgb (g/dL) 9.90 (8.50, 11.40) 9.80 (8.40, 11.50) 0.9 10.00 (8.70, 11.70) 9.60 (8.25, 11.20) 0.027
Baseline ferritin (ng/mL) 346 (140, 1,250) 379 (125, 1,219) >0.9 293 (137, 832) 393 (128, 1,254) 0.3
Heavy chain

0.001

<0.001
    IgG 72 (58%) 109 (44%)
72 (59%) 107 (43%)
    IgA 26 (21%) 43 (17%)
26 (21%) 42 (17%)
    IgM 1 (0.7%) 1 (0.4%)
1 (0.8%) 1 (0.4%)
    IgD 2 (1.5%) 0 (0%)
2 (1.6%) 0 (0%)
    None 23 (19%) 95 (38%)
22 (18%) 98 (40%)
True EMD 20 (16%) 39 (16%) >0.9 26 (21%) 34 (14%) 0.067
≥1 HRCA 81 (70%) 151 (65%) 0.3 81 (70%) 152 (65%) 0.4
Prior LOTs 6.00 (5.00, 8.00) 6.00 (4.00, 8.00) 0.4 6.00 (4.00, 8.00) 6.00 (4.00, 8.00) >0.9
Prior ASCT 66 (54%) 141 (57%) 0.6 61 (50%) 144 (58%) 0.12
Prior CAR-T 53 (43%) 90 (36%) 0.2 50 (41%) 94 (38%) 0.6
Triple refractory 112 (91%) 224 (91%) 0.9 112 (91%) 226 (91%) >0.9
Penta refractory 61 (49%) 123 (49%) >0.9 61 (50%) 126 (51%) 0.8
Trial eligible 24 (19%) 55 (22%) 0.5 26 (21%) 53 (21%) >0.9
Prior BCMA-directed therapy 67 (54%) 133 (54%) >0.9 60 (49%) 138 (56%) 0.2
Received full dose 114 (93%) 241 (97%) 0.059 114 (93%) 241 (97%) 0.045
1 Median (Q1, Q3); n (%)
2 Design-based KruskalWallis test; Pearson’s X^2: Rao & Scott adjustment
3 Median (Min to Max); n (%); Median (Q1, Q3)
4 Wilcoxon rank sum test; Pearson’s Chi-squared test; Fisher’s exact test

5 STACKED BAR GRAPH

5.1 FIGURE 1A: Toxicity, IPTW

data1 = IPTW.data %>%
  transmute(
    BsAb = as.factor(bs_ab_name),
    CRS.any = ifelse( max_crs_grade_0_5 != 0, 1, 0),
    ICANS.any = ifelse( max_icans_grade != 0,1,0),
    
    CRS.two = ifelse( max_crs_grade_0_5 == 2 | max_crs_grade_0_5 == 3 | max_crs_grade_0_5 == 4, 1, 0),
    ICANS.two = ifelse( max_icans_grade ==2 | max_icans_grade == 3 | max_icans_grade == 4,1,0),
    
    CRS.three = ifelse( max_crs_grade_0_5 == 3 | max_crs_grade_0_5 == 4, 1, 0),
    ICANS.three = ifelse( max_icans_grade == 3 | max_icans_grade == 4,1,0),
    weights
  )


weighted_table <- survey::svydesign(
    id = ~0,
    weights = ~weights,
    data = data1,
    fpc = NULL
  )


combined_data <- rbind(
  IPTW.data %>%
    count(bs_ab_name, max_crs_grade_0_5, wt=weights) %>%
    group_by(bs_ab_name) %>%
    mutate(
      bs_ab_name = as.factor(recode(bs_ab_name, "tec"="Tec","elra"="Elra")),
      total = sum(n[max_crs_grade_0_5 != 0]),
      pct = prop.table(n) * 100,
      total_pct = sum(pct[max_crs_grade_0_5 != 0]),
      Grade = factor(max_crs_grade_0_5, levels = c(0, 4, 3, 2, 1)),
      Category = "CRS"
    ),
  IPTW.data %>%
    count(bs_ab_name, max_icans_grade, wt=weights) %>%
    group_by(bs_ab_name) %>%
    mutate(
      bs_ab_name = as.factor(recode(bs_ab_name, "tec"="Tec","elra"="Elra")),
      total = sum(n[max_icans_grade != 0]),
      pct = prop.table(n) * 100,
      total_pct = sum(pct[max_icans_grade != 0]),
      Grade = factor(max_icans_grade, levels = c(0, 4, 3, 2, 1)),
      Category = "ICANS"
    )
) %>%
  mutate(
    Category = factor(Category, levels = c("CRS","ICANS"))
  )


# Any grade toxicity
CRS.pval <- as.numeric(svychisq(~CRS.any+BsAb, weighted_table)$p.value)
ICANS.pval <- as.numeric(svychisq(~ICANS.any+BsAb, weighted_table)$p.value)

f_labels <- data.frame(Category = c("CRS","ICANS"), label = c(CRS.pval, ICANS.pval) ) %>%
  mutate(
    Category = factor(Category, levels = c("CRS","ICANS"))
  )

# G2+ toxicity
CRS.two.pval <- as.numeric(svychisq(~CRS.two+BsAb, weighted_table)$p.value)
ICANS.two.pval <- as.numeric(svychisq(~ICANS.two+BsAb, weighted_table)$p.value)

f_labels_two <- data.frame(Category = c("CRS","ICANS"), label = c(CRS.two.pval, ICANS.two.pval) ) %>%
  mutate(
    Category = factor(Category, levels = c("CRS","ICANS"))
  )

# G3+ toxicity
CRS.three.pval <- as.numeric(svychisq(~CRS.three+BsAb, weighted_table)$p.value)
ICANS.three.pval <- as.numeric(svychisq(~ICANS.three+BsAb, weighted_table)$p.value)

f_labels_three <- data.frame(Category = c("CRS","ICANS"), label = c(CRS.three.pval, ICANS.three.pval) ) %>%
  mutate(
    Category = factor(Category, levels = c("CRS","ICANS"))
  )




toxplot_IPTW <- ggplot(combined_data) +
  aes(x = bs_ab_name, y = pct, fill = Grade) +
  geom_bar(stat = "identity", position = "stack", data = . %>% filter(Grade != 0)) +
    geom_text(aes(label = paste0(signif(pct,2), "% (", round(n,0) ,")")), 
            data = . %>% filter(Grade != 0 & Grade != 4),
            position = position_stack(vjust = 0.5)) + 
  geom_text( aes(label = paste0(signif(total_pct,2), "% (", signif(total,2),")"), x = bs_ab_name, y = total_pct + 6, vjust = 0), 
             data = . %>% filter(Grade == 1),
             color = "#104a8e",
             fontface = "bold") +
  facet_wrap(~Category, scales = "free_y") +
  ylab("Proportion of patients (% [n])") +
  xlab(NULL) +
  theme_classic()+
  scale_fill_brewer(palette = "Pastel1") +
  #guides(fill = guide_legend(reverse = TRUE)) +
  scale_y_continuous(breaks = seq(0, 100, by = 20), limits = c(0, 110)) +
  theme(
    plot.title = element_text(size = 14),
    axis.title = element_text(size = 14),
    axis.text = element_text(size = 14),
    legend.title = element_text(size = 14),
    legend.text = element_text(size = 14),
    strip.text = element_text(size = 14),
    axis.text.x = element_text(angle = 45, hjust = 1),
    #legend.position = "bottom"
    legend.position = "right"
  ) + 
  coord_cartesian(ylim = c(0, 115)) +
  geom_text( aes(label = paste0("Grade 1+, p=", signif(label,2)) ), x = 1.5, y = 100+14, vjust = 0, inherit.aes = FALSE, data = f_labels ) + 
  geom_text( aes(label = paste0("Grade 2+, p=", signif(label,2)) ), x = 1.5, y = 100+7, vjust = 0, inherit.aes = FALSE, data = f_labels_two ) + 
  geom_text( aes(label = paste0("Grade 3+, p=", signif(label,2)) ), x = 1.5, y = 100, vjust = 0, inherit.aes = FALSE, data = f_labels_three ) + 
  ggtitle("A) IPTW")

toxplot_IPTW

5.2 FIGURE 1B: Toxicity, unweighted

combined_data <- rbind(
  combined_db %>%
    count(bs_ab_name, max_crs_grade_0_5) %>%
    group_by(bs_ab_name) %>%
    mutate(
      total = sum(n[max_crs_grade_0_5 != 0]),
      pct = prop.table(n) * 100,
      total_pct = sum(pct[max_crs_grade_0_5 != 0]),
      Grade = factor(max_crs_grade_0_5, levels = c(0, 4, 3, 2, 1)),
      bs_ab_name = as.factor(recode(bs_ab_name, "tec"="Tec","elra"="Elra")),
      Category = "CRS"
    ),
  combined_db %>%
    count(bs_ab_name, max_icans_grade) %>%
    group_by(bs_ab_name) %>%
    mutate(
      total = sum(n[max_icans_grade != 0]),
      pct = prop.table(n) * 100,
      total_pct = sum(pct[max_icans_grade != 0]),
      Grade = factor(max_icans_grade, levels = c(0, 4, 3, 2, 1)),
      bs_ab_name = as.factor(recode(bs_ab_name, "tec"="Tec","elra"="Elra")),
      Category = "ICANS"
    )
) %>%
  mutate(
    Category = factor(Category, levels = c("CRS","ICANS"))
    
  )




## Any-grade toxicity
CRS.d.test <- data.frame(
  results = ifelse(combined_db$max_crs_grade_0_5 !=0,1,0) ,
  test = combined_db$bs_ab_name
)
CRS.pval <- chisq.test(CRS.d.test$results,CRS.d.test$test, correct = FALSE)$p.value

ICANS.d.test <- data.frame(
  results = ifelse(combined_db$max_icans_grade !=0,1,0) ,
  test = combined_db$bs_ab_name
)
ICANS.pval <- chisq.test(ICANS.d.test$results,ICANS.d.test$test, correct = FALSE)$p.value


f_labels <- data.frame(Category = c("CRS","ICANS"), label = c(CRS.pval, ICANS.pval) ) %>%
  mutate(
    Category = factor(Category, levels = c("CRS","ICANS") )
  )


## G2+ toxicity
CRS.two.d.test <- data.frame(
  results = ifelse(combined_db$max_crs_grade_0_5 == 2 | combined_db$max_crs_grade_0_5 == 3 | combined_db$max_crs_grade_0_5 == 4,1,0),
  test = combined_db$bs_ab_name
)
CRS.two.pval <- chisq.test(CRS.two.d.test$results,CRS.two.d.test$test, correct = FALSE)$p.value

ICANS.two.d.test <- data.frame(
  results = ifelse(combined_db$max_icans_grade == 2 | combined_db$max_icans_grade == 3 | combined_db$max_icans_grade == 4,1,0),
  test = combined_db$bs_ab_name
)
ICANS.two.pval <- chisq.test(ICANS.two.d.test$results,ICANS.two.d.test$test, correct = FALSE)$p.value

f_labels_two <- data.frame(Category = c("CRS","ICANS"), label = c(CRS.two.pval, ICANS.two.pval) ) %>%
  mutate(
    Category = factor(Category, levels = c("CRS","ICANS") )
  )

## G3+ toxicity
CRS.three.d.test <- data.frame(
  results = ifelse( combined_db$max_crs_grade_0_5 == 3 | combined_db$max_crs_grade_0_5 == 4,1,0),
  test = combined_db$bs_ab_name
)
CRS.three.pval <- fisher.test(CRS.three.d.test$results,CRS.three.d.test$test)$p.value

ICANS.three.d.test <- data.frame(
  results = ifelse(combined_db$max_icans_grade == 3 | combined_db$max_icans_grade == 4,1,0),
  test = combined_db$bs_ab_name
)
ICANS.three.pval <- chisq.test(ICANS.three.d.test$results,ICANS.three.d.test$test, correct = FALSE)$p.value

f_labels_three <- data.frame(Category = c("CRS","ICANS"), label = c(CRS.three.pval, ICANS.three.pval) ) %>%
  mutate(
    Category = factor(Category, levels = c("CRS","ICANS") )
  )



toxplot_unweighted <- ggplot(combined_data) +
  aes(x = bs_ab_name, y = pct, fill = Grade) +
  geom_bar(stat = "identity", position = "stack", data = . %>% filter(Grade != 0)) +
    geom_text(aes(label = paste0(signif(pct,2), "% (",n,")")), 
            data = . %>% filter(Grade != 0 & Grade != 4),
            position = position_stack(vjust = 0.5)) + 
  geom_text( aes(label = paste0(signif(total_pct,2), "% (",total,")"), x = bs_ab_name, y = total_pct + 6, vjust = 0), 
             data = . %>% filter(Grade == 1),
             color = "#104a8e",
             fontface = "bold") +
  facet_wrap(~Category, scales = "free_y") +
  ylab("Proportion of patients (% [n])") +
  xlab(NULL) +
  theme_classic()+
  scale_y_continuous(breaks = seq(0, 100, by = 20), limits = c(0, 100)) +
  theme(
    plot.title = element_text(size = 14),
    axis.title = element_text(size = 14),
    axis.text = element_text(size = 14),
    legend.title = element_text(size = 14),
    legend.text = element_text(size = 14),
    strip.text = element_text(size = 14),
    axis.text.x = element_text(angle = 45, hjust = 1),
    legend.position = "none"
  ) + 
  coord_cartesian(ylim = c(0, 114)) +
  scale_fill_brewer(palette = "Pastel1", direction = 1) +
  geom_text( aes(label = paste0("Grade 1+, p=", signif(label,2)) ), x = 1.5, y = 100+14, vjust = 0, inherit.aes = FALSE, data = f_labels ) + 
  geom_text( aes(label = paste0("Grade 2+, p=", signif(label,2)) ), x = 1.5, y = 100+7, vjust = 0, inherit.aes = FALSE, data = f_labels_two ) + 
  geom_text( aes(label = paste0("Grade 3+, p=", signif(label,2)) ), x = 1.5, y = 100, vjust = 0, inherit.aes = FALSE, data = f_labels_three ) + 
  ggtitle("B) Unweighted")

toxplot_unweighted

5.3 ***FIGURE 1: IPTW and unweighted toicity

FIGURE1 <- toxplot_IPTW / toxplot_unweighted

ggsave("svg/FIGURE1.svg", FIGURE1, device = "svg",
       width = 7,
       height = 11,
       units = c("in"))

ggsave("tiff/FIGURE1.tiff", FIGURE1, device = "tiff",
       width = 7,
       height = 11,
       dpi = 600,
       units = c("in"),
       compression = "lzw")

ggsave("jpeg/FIGURE1.jpeg", FIGURE1, device = "jpeg",
       width = 7,
       height = 11,
       dpi = 600,
       units = "in")

knitr::include_graphics("jpeg/FIGURE1.jpeg")

5.4 FIGURE 2A: Response, IPTW

data1 = IPTW.data %>%
  filter(!is.na(best_response)) %>%
  transmute(
    BsAb = as.factor(bs_ab_name),
    best_response = factor(best_response, levels =  c("sCR", "CR","VGPR","PR","MR","SD","PD")),
    ORR = ifelse(best_response == "sCR" | best_response == "CR" | best_response == "VGPR" | best_response == "PR",1,0),
    VGPR = ifelse(best_response == "sCR" | best_response == "CR" | best_response == "VGPR",1,0),
    CR = ifelse(best_response == "sCR" | best_response == "CR",1,0),
    weights
  )

weighted_table <- survey::svydesign(
    id = ~0,
    weights = ~weights,
    data = data1,
    fpc = NULL
  )



data1 <-   IPTW.data %>%
  filter(!is.na(best_response)) %>%
  count(bs_ab_name, best_response, wt=weights) %>%
  group_by(bs_ab_name) %>%
  mutate(
    total = sum(n[best_response == "sCR" | best_response == "CR" | best_response == "VGPR" | best_response == "PR"]),
    pct = prop.table(n) * 100,
    total_pct = sum(pct[best_response == "sCR" | best_response == "CR" | best_response == "VGPR" | best_response == "PR"]),
    best_response = factor(best_response, levels =  c("sCR", "CR","VGPR","PR","MR", "SD","PD")),
    bs_ab_name = as.factor(bs_ab_name),
    Category = "Response rate"
  )



Response.pval <- as.numeric(svychisq(~ORR+BsAb, weighted_table)$p.value)
CR.pval <- as.numeric(svychisq(~CR+BsAb, weighted_table)$p.value)
VGPR.pval <- as.numeric(svychisq(~VGPR+BsAb, weighted_table)$p.value)

f_labels <- data.frame(Category = c("Response rate"), label = c(Response.pval) ) %>%
  mutate(
    Category = factor(Category, levels = c("Response rate"))
  )

f_labels_CR <- data.frame(Category = c("Response rate"), label = c(CR.pval) ) %>%
  mutate(
    Category = factor(Category, levels = c("Response rate"))
  )

f_labels_VGPR <- data.frame(Category = c("Response rate"), label = c(VGPR.pval) ) %>%
  mutate(
    Category = factor(Category, levels = c("Response rate"))
  )

responseplot_IPTW <- ggplot(data1) +
  aes(x = bs_ab_name, y = pct, fill = best_response) +
  geom_bar(stat = "identity", position = "stack", data = . %>% filter(best_response == "sCR" | best_response == "CR" | best_response == "VGPR" | best_response == "PR")) +
    geom_text(aes(label = paste0(signif(pct,2), "% (", round(n,0) ,")")), 
            data = . %>% filter(best_response == "sCR" | best_response == "CR" | best_response == "VGPR"| best_response == "PR"),
            position = position_stack(vjust = 0.5)) + 
  geom_text( aes(label = paste0(signif(total_pct,2), "% (", signif(total,2),")"), x = bs_ab_name, y = total_pct + 7, vjust = 0), 
             data = . %>% filter(best_response == "PR"),
             color = "#104a8e",
             fontface = "bold") +
  facet_wrap(~Category, scales = "free_y") +
  ylab("Proportion of patients (% [n])") +
  xlab(NULL) +
  theme_classic()+
  scale_fill_brewer(
    palette = "Pastel1",
    direction = -1,
    breaks = rev(levels(data1$best_response))
  ) +
  guides(fill = guide_legend(reverse = TRUE)) +
  scale_y_continuous(breaks = seq(0, 100, by = 20), limits = c(0, 110)) +
  theme(
    plot.title = element_text(size = 14),
    axis.title = element_text(size = 14),
    axis.text = element_text(size = 14),
    legend.title = element_text(size = 14),
    legend.text = element_text(size = 14),
    strip.text = element_text(size = 14),
    axis.text.x = element_text(angle = 45, hjust = 1),
    #legend.position = "bottom"
    legend.position = "right"
  ) + 
  coord_cartesian(ylim = c(0, 114)) +
  geom_text( aes(label = paste0("ORR, p=", signif(Response.pval,2)) ), x = 1.5, y = 100+14, vjust = 0, inherit.aes = FALSE, data = f_labels  ) + 
  geom_text( aes(label = paste0("≥VGPR, p=", signif(VGPR.pval,2)) ), x = 1.5, y = 100+7, vjust = 0, inherit.aes = FALSE, data = f_labels_VGPR  ) + 
  geom_text( aes(label = paste0("≥CR, p=", signif(CR.pval,2)) ), x = 1.5, y = 100, vjust = 0, inherit.aes = FALSE, data = f_labels_CR  ) + 
  labs(fill = "Best response")+ 
  ggtitle("A) IPTW")

responseplot_IPTW

5.5 FIGURE 2B: Response, unweighted

data1 <- combined_db %>%
  filter(!is.na(best_response)) %>%
  transmute(
    best_response = factor(
      best_response, 
      levels =  c("sCR", "CR","VGPR", "PR","MR","SD","PD")),
    ORR = ifelse(best_response == "sCR" | best_response == "CR" | best_response == "VGPR" |best_response == "PR",1,0),
    bs_ab_name = as.factor(recode(bs_ab_name, "tec"="Tec","elra"="Elra")),
    Category = "Response rate"
  )

combined_data <- combined_db %>%
  filter(!is.na(best_response)) %>%
  count(bs_ab_name, best_response) %>%
  group_by(bs_ab_name) %>%
  transmute(
    bs_ab_name = as.factor(recode(bs_ab_name, "tec"="Tec","elra"="Elra")),
    best_response = factor(
      best_response, 
      levels =  c("sCR", "CR","VGPR", "PR","MR","SD","PD")),
    total = sum(n[best_response == "sCR" | best_response == "CR" | best_response == "VGPR" | best_response == "PR"]),
    pct = prop.table(n) * 100,
    n,
    total_pct = sum(pct[best_response == "sCR" | best_response == "CR" |best_response == "VGPR" | best_response == "PR"]),
    Category = "Response rate"
  )

## ORR
Response.d.test <- data.frame(
  results = ifelse(data1$ORR,1,0),
  test = data1$bs_ab_name
)
Response.pval <- chisq.test(Response.d.test$results,Response.d.test$test, correct = FALSE)$p.value

## VGPR
VGPR.d.test <- data.frame(
  results = ifelse(data1$best_response == "CR" | data1$best_response == "VGPR",1,0),
  test = data1$bs_ab_name
)
VGPR.pval <- chisq.test(VGPR.d.test$results,VGPR.d.test$test, correct = FALSE)$p.value

## CR
CR.d.test <- data.frame(
  results = ifelse(data1$best_response == "CR",1,0),
  test = data1$bs_ab_name
)
CR.pval <- chisq.test(CR.d.test$results,CR.d.test$test, correct = FALSE)$p.value


f_labels <- data.frame(Category = c("Response rate"), label = c(Response.pval) ) %>%
  mutate(
    Category = factor(Category, levels = c("Response rate"))
  )

f_labels_CR <- data.frame(Category = c("Response rate"), label = c(CR.pval) ) %>%
  mutate(
    Category = factor(Category, levels = c("Response rate"))
  )

f_labels_VGPR <- data.frame(Category = c("Response rate"), label = c(VGPR.pval) ) %>%
  mutate(
    Category = factor(Category, levels = c("Response rate"))
  )


responseplot_unweighted <- ggplot(combined_data) +
  aes(x = bs_ab_name, y = pct, fill = best_response) +
  geom_bar(stat = "identity", position = "stack", data = . %>% filter(best_response == "sCR" | best_response == "CR" |best_response == "VGPR" | best_response == "PR")) +
  geom_text(aes(label = paste0(signif(pct,2), "% (",n,")")), 
             data = . %>% filter(best_response == "sCR" | best_response == "CR" | best_response == "VGPR" | best_response == "PR"),
             position = position_stack(vjust = 0.5)) + 
  geom_text( aes(label = paste0(signif(total_pct,2), "% (",total,")"), x = bs_ab_name, y = total_pct + 7, vjust = 0), 
             data = . %>% filter(best_response == "PR"),
             color = "#104a8e",
             fontface = "bold") +
  facet_wrap(~Category, scales = "free_y") +
  ylab("Proportion of patients (% [n])") +
  xlab(NULL) +
  theme_classic()+
  scale_y_continuous(breaks = seq(0, 100, by = 20), limits = c(0, 110)) +
  theme(
    plot.title = element_text(size = 14),
    axis.title = element_text(size = 14),
    axis.text = element_text(size = 14),
    legend.title = element_text(size = 14),
    legend.text = element_text(size = 14),
    strip.text = element_text(size = 14),
    axis.text.x = element_text(angle = 45, hjust = 1),
    #legend.position = "bottom"
    legend.position = "none"
  ) + 
  coord_cartesian(ylim = c(0, 114)) +
  scale_fill_brewer(palette = "Pastel1", direction = -1) +
  geom_text( aes(label = paste0("ORR, p=", signif(Response.pval,2)) ), x = 1.5, y = 100+14, vjust = 0, inherit.aes = FALSE, data = f_labels) + 
  geom_text( aes(label = paste0("≥VGPR, p=", signif(VGPR.pval,2)) ), x = 1.5, y = 100+7, vjust = 0, inherit.aes = FALSE, data = f_labels_VGPR  ) + 
  geom_text( aes(label = paste0("≥CR, p=", signif(CR.pval,2)) ), x = 1.5, y = 100, vjust = 0, inherit.aes = FALSE, data = f_labels_CR  ) + 
  labs(fill = "Best response")+ 
  ggtitle("B) Unweighted")

responseplot_unweighted

5.6 ***FIGURE 2: IPTW and unweighted response

FIGURE2 <- responseplot_IPTW / responseplot_unweighted

ggsave("svg/FIGURE2.svg", FIGURE2, device = "svg",
       width = 7,
       height = 11,
       units = c("in"))

ggsave("tiff/FIGURE2.tiff", FIGURE2, device = "tiff",
       width = 7,
       height = 11,
       dpi = 600,
       units = c("in"),
       compression = "lzw")

ggsave("jpeg/FIGURE2.jpeg", FIGURE2, device = "jpeg",
       width = 7,
       height = 11,
       dpi = 600,
       units = "in")

knitr::include_graphics("jpeg/FIGURE2.jpeg")

6 SURVIVAL ANALYSIS

6.1 Length of follow-up (Reverse KM, months)

data1 <- IPTW.data %>%
  transmute(
    bs_ab_name,
    Death = ifelse(death_yes_no,1,0),
    Days.to.death.or.DLC = as.numeric(date_of_last_contact - date_of_bs_ab_step_up_d1),
    Reverse_death = ifelse(death_yes_no == 1, 0,1),
    weights
  )

reverse_km_OS <- survfit(Surv(Days.to.death.or.DLC/30, Reverse_death) ~ 1, weights = weights, data1)
reverse_km_OS
## Call: survfit(formula = Surv(Days.to.death.or.DLC/30, Reverse_death) ~ 
##     1, data = data1, weights = weights)
## 
##    15 observations deleted due to missingness 
##      records   n events median 0.95LCL 0.95UCL
## [1,]     356 356    215   13.4    12.5    14.9
reverse_km_OS_elra <- survfit(Surv(Days.to.death.or.DLC/30, Reverse_death) ~ 1, weights = weights, data1 %>% filter(bs_ab_name == "elra"))
reverse_km_OS_elra
## Call: survfit(formula = Surv(Days.to.death.or.DLC/30, Reverse_death) ~ 
##     1, data = data1 %>% filter(bs_ab_name == "elra"), weights = weights)
## 
##    6 observations deleted due to missingness 
##      records   n events median 0.95LCL 0.95UCL
## [1,]     117 117   74.5   7.47    6.03       9
reverse_km_OS_tec <- survfit(Surv(Days.to.death.or.DLC/30, Reverse_death) ~ 1, weights = weights, data1 %>% filter(bs_ab_name == "tec"))
reverse_km_OS_tec
## Call: survfit(formula = Surv(Days.to.death.or.DLC/30, Reverse_death) ~ 
##     1, data = data1 %>% filter(bs_ab_name == "tec"), weights = weights)
## 
##    9 observations deleted due to missingness 
##      records   n events median 0.95LCL 0.95UCL
## [1,]     239 239    140   17.1    14.9    20.1

6.2 Stratified by BsAb

data1 <- IPTW.data %>%
  transmute(
    site,
    patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate,
    bs_ab_name,
    Relapse.or.death = ifelse(progression_yes_no | death_yes_no,1,0),
    Death = ifelse(death_yes_no,1,0),
    date_of_bs_ab_step_up_d1,
    date_of_pd,
    date_of_last_contact,
    date_progression_or_death = coalesce(date_of_pd, date_of_last_contact),
    Days.to.death.or.DLC = as.numeric(date_of_last_contact - date_of_bs_ab_step_up_d1),
    Days.to.death.or.relapse.or.DLC = as.numeric(date_progression_or_death  - date_of_bs_ab_step_up_d1),
    trial_eligible = ifelse(eligible_for_corresponding_bs_ab_trial_on_first_dose_use_trial_elig_sheet_to_help_you ==1,1,0),
    weights
  )

data_dor <- IPTW.data %>%
  filter(best_response %in% c("sCR","CR","VGPR","PR")) %>%
  transmute(
    weights,
    bs_ab_name,
    day_30_response,
    day_90_response,
    best_response,
    Relapse.or.death = ifelse(progression_yes_no | death_yes_no,1,0),
    Death = ifelse(death_yes_no,1,0),
    date_of_pd,
    date_of_last_contact,
    date_of_bs_ab_step_up_d1,
    date_progression_or_death = coalesce(date_of_pd, date_of_last_contact),
    day_30_response,
    day_90_response,
    best_response,
    time_to_best_response = case_when(
      best_response == day_30_response ~ 30,
      best_response == day_90_response ~ 90,
      .default = NA
    ),
    trial_eligible = ifelse(eligible_for_corresponding_bs_ab_trial_on_first_dose_use_trial_elig_sheet_to_help_you ==1,1,0)
  ) %>%
  filter(!is.na(time_to_best_response)) %>%
  mutate(
    Days.to.death.or.DLC = as.numeric(date_of_last_contact - date_of_bs_ab_step_up_d1) - time_to_best_response,
    Days.to.death.or.relapse.or.DLC = as.numeric(date_progression_or_death  - date_of_bs_ab_step_up_d1)- time_to_best_response,
  ) %>%
  filter(Days.to.death.or.DLC >0 & Days.to.death.or.relapse.or.DLC >0)

km_OS <- survfit(Surv(Days.to.death.or.DLC/30, Death) ~ bs_ab_name, weights = weights, data = data1)
km_OS
## Call: survfit(formula = Surv(Days.to.death.or.DLC/30, Death) ~ bs_ab_name, 
##     data = data1, weights = weights)
## 
##    15 observations deleted due to missingness 
##                 records   n events median 0.95LCL 0.95UCL
## bs_ab_name=elra     117 117   42.8     NA     8.6      NA
## bs_ab_name=tec      239 239   98.6   21.5    16.4      NA
med_time <- surv_median(km_OS)

cox_model <- coxph(Surv(Days.to.death.or.DLC/30, Death) ~ ifelse(bs_ab_name == "elra",1,0), weights = weights, data = data1)
cox_summary <- summary(cox_model)

HR <- round(cox_summary$coefficients[1, "exp(coef)"], 2)
CI_lower <- round(cox_summary$conf.int[1, "lower .95"], 2)
CI_upper <- round(cox_summary$conf.int[1, "upper .95"], 2)
pval <- signif(cox_summary$coefficients[1, "Pr(>|z|)"], 2)

## 1-year timepoint
timepoint1 <- summary(km_OS,times=c(3))
timepoint2 <- summary(km_OS,times=c(6))
timepoint3 <- summary(km_OS,times=c(9))
timepoint4 <- summary(km_OS,times=c(12))
timepoint5 <- summary(km_OS,times=c(15))

OS <- ggsurvplot(km_OS,
           pval=FALSE,
           conf.int = FALSE,
           risk.table = TRUE, # Add risk table
           fontsize = 4,
           risk.table.col = "strata", # Change risk table color by groups
           tables.height = 0.3,
           tables.theme = theme_cleantable() +
             theme(
               axis.text.y = element_text(size = 9),
               plot.title = element_text(size = 10)
               ),
           xlab="Time (months)", ylab = "Overall survival",
           xlim = c(0,18), break.x.by = c(3),
           ylim = c(0,1), break.y.by = c(0.25),
           legend = "none",
           surv.scale="percent",
           font.main = c(10, "plain", "black"),
           font.x = c(10, "plain", "black"),
           font.y = c(10, "plain", "black"),
           font.caption = c(10, "plain", "black"),
           font.tickslab = c(10, "plain", "black"),
           legend.labs = c("Elra", "Tec"),
           title = "A)",
           size = 0.5
) 

OS$plot <- OS$plot + 
  geom_vline(xintercept = c(3, 6, 9, 12, 15), linetype = "dotted", color = "gray50") +
  annotate("text", x = 3, y = timepoint1$surv[1] - 0.1, label = paste0(signif(timepoint1$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 3, y = timepoint1$surv[2] + 0.03, label = paste0(signif(timepoint1$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 6, y = timepoint2$surv[1] - 0.08, label = paste0(signif(timepoint2$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 6, y = timepoint2$surv[2] + 0.03, label = paste0(signif(timepoint2$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 9, y = timepoint3$surv[1] - 0.08, label = paste0(signif(timepoint3$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 9, y = timepoint3$surv[2] + 0.03, label = paste0(signif(timepoint3$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 12, y = timepoint4$surv[1] - 0.08, label = paste0(signif(timepoint4$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 12, y = timepoint4$surv[2] + 0.03, label = paste0(signif(timepoint4$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 15, y = timepoint5$surv[1] - 0.08, label = paste0(signif(timepoint5$surv[1]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 15, y = timepoint5$surv[2] + 0.03, label = paste0(signif(timepoint5$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = -0.7, y = 0,
           label = paste0("Elra vs tec\n",
                          " Median: ", signif(med_time$median[1], 3), " vs ", signif(med_time$median[2], 3), " months\n",
                          " HR: ", sprintf("%.2f", HR), " (95% CI ", sprintf("%.2f", CI_lower), "-", sprintf("%.2f", CI_upper),", p=",pval, ")"),
           color = "black", hjust = 0, vjust = 0, size = 3)

OS

km_PFS <- survfit(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ bs_ab_name, weights = weights, data = data1)
km_PFS
## Call: survfit(formula = Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ 
##     bs_ab_name, data = data1, weights = weights)
## 
##    12 observations deleted due to missingness 
##                 records   n events median 0.95LCL 0.95UCL
## bs_ab_name=elra     118 118   68.1   3.77    2.57    7.17
## bs_ab_name=tec      241 241  151.1   8.50    7.10   11.30
med_time <- surv_median(km_PFS)

cox_model <- coxph(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ ifelse(bs_ab_name == "elra",1,0), weights = weights, data = data1)
cox_summary <- summary(cox_model)

HR <- round(cox_summary$coefficients[1, "exp(coef)"], 2)
CI_lower <- round(cox_summary$conf.int[1, "lower .95"], 2)
CI_upper <- round(cox_summary$conf.int[1, "upper .95"], 2)
pval <- signif(cox_summary$coefficients[1, "Pr(>|z|)"], 2)

## Timepoints
timepoint1 <- summary(km_PFS,times=c(3))
timepoint2 <- summary(km_PFS,times=c(6))
timepoint3 <- summary(km_PFS,times=c(9))
timepoint4 <- summary(km_PFS,times=c(12))
timepoint5 <- summary(km_PFS,times=c(15))


PFS <- ggsurvplot(km_PFS,
           pval=FALSE,
           conf.int = FALSE,
           risk.table = TRUE, # Add risk table
           fontsize = 4,
           risk.table.col = "strata", # Change risk table color by groups
           tables.height = 0.3,
           tables.theme = theme_cleantable() +
             theme(
               axis.text.y = element_text(size = 9),
               plot.title = element_text(size = 10)
               ),
           xlab="Time (months)", ylab = "Progression-free survival",
           xlim = c(0,18), break.x.by = c(3),
           ylim = c(0,1), break.y.by = c(0.25),
           legend = "none",
           surv.scale="percent",
           font.main = c(10, "plain", "black"),
           font.x = c(10, "plain", "black"),
           font.y = c(10, "plain", "black"),
           font.caption = c(10, "plain", "black"),
           font.tickslab = c(10, "plain", "black"),
           legend.labs = c("Elra", "Tec"),
           title = "B)",
           size = 0.5
)

PFS$plot <- PFS$plot + 
  geom_vline(xintercept = c(3, 6, 9, 12, 15), linetype = "dotted", color = "gray50") +
  annotate("text", x = 3, y = timepoint1$surv[1] - 0.08, label = paste0(signif(timepoint1$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 3, y = timepoint1$surv[2] + 0.03, label = paste0(signif(timepoint1$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 6, y = timepoint2$surv[1] - 0.08, label = paste0(signif(timepoint2$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 6, y = timepoint2$surv[2] + 0.03, label = paste0(signif(timepoint2$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 9, y = timepoint3$surv[1] - 0.08, label = paste0(signif(timepoint3$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 9, y = timepoint3$surv[2] + 0.03, label = paste0(signif(timepoint3$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 12, y = timepoint4$surv[1] - 0.08, label = paste0(signif(timepoint4$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 12, y = timepoint4$surv[2] + 0.03, label = paste0(signif(timepoint4$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 15, y = timepoint5$surv[1] - 0.08, label = paste0(signif(timepoint5$surv[1]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 15, y = timepoint5$surv[2] + 0.03, label = paste0(signif(timepoint5$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = -0.7, y = 0,
           label = paste0("Elra vs tec\n",
                          " Median: ", signif(med_time$median[1], 3), " vs ", signif(med_time$median[2], 3), " months\n",
                          " HR: ", sprintf("%.2f", HR), " (95% CI ", sprintf("%.2f", CI_lower), "-", sprintf("%.2f", CI_upper),", p=",pval, ")"),
           color = "black", hjust = 0, vjust = 0, size = 3)

PFS

km_DOR <- survfit(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ bs_ab_name, weights = weights, data = data_dor)
km_DOR
## Call: survfit(formula = Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ 
##     bs_ab_name, data = data_dor, weights = weights)
## 
##    3 observations deleted due to missingness 
##                 records     n events median 0.95LCL 0.95UCL
## bs_ab_name=elra      62  60.2   20.6   11.7    6.13      NA
## bs_ab_name=tec      118 120.7   56.1   13.7    9.27    22.3
med_time <- surv_median(km_DOR)

cox_model <- coxph(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ ifelse(bs_ab_name == "elra",1,0), weights = weights, data = data_dor)
cox_summary <- summary(cox_model)

HR <- round(cox_summary$coefficients[1, "exp(coef)"], 2)
CI_lower <- round(cox_summary$conf.int[1, "lower .95"], 2)
CI_upper <- round(cox_summary$conf.int[1, "upper .95"], 2)
pval <- signif(cox_summary$coefficients[1, "Pr(>|z|)"], 2)

## Specific timepoints
timepoint1 <- summary(km_DOR,times=c(3))
timepoint2 <- summary(km_DOR,times=c(6))
timepoint3 <- summary(km_DOR,times=c(9))
timepoint4 <- summary(km_DOR,times=c(12))
timepoint5 <- summary(km_DOR,times=c(15))


DOR <- ggsurvplot(km_DOR,
           pval=FALSE,
           conf.int = FALSE,
           risk.table = TRUE, # Add risk table
           fontsize = 4,
           risk.table.col = "strata", # Change risk table color by groups
           tables.height = 0.3,
           #risk.table.fontsize = 4,
           tables.theme = theme_cleantable() +
             theme(
               axis.text.y = element_text(size = 9),
               plot.title = element_text(size = 10)
               ),
           xlab="Time (months)", ylab = "Duration of response",
           xlim = c(0,18), break.x.by = c(3),
           ylim = c(0,1), break.y.by = c(0.25),
           legend = "none",
           surv.scale="percent",
           font.main = c(10, "plain", "black"),
           font.x = c(10, "plain", "black"),
           font.y = c(10, "plain", "black"),
           font.caption = c(10, "plain", "black"),
           font.tickslab = c(10, "plain", "black"),
           legend.labs = c("Elra", "Tec"),
           title = "C)",
           size = 0.5
)

DOR$plot <- DOR$plot + 
  geom_vline(xintercept = c(3, 6, 9, 12, 15), linetype = "dotted", color = "gray50") +
  annotate("text", x = 3, y = timepoint1$surv[1] -0.08, label = paste0(signif(timepoint1$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 3, y = timepoint1$surv[2] + 0.03, label = paste0(signif(timepoint1$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 6, y = timepoint2$surv[1] - 0.08, label = paste0(signif(timepoint2$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 6, y = timepoint2$surv[2] + 0.03, label = paste0(signif(timepoint2$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 9, y = timepoint3$surv[1] - 0.08, label = paste0(signif(timepoint3$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 9, y = timepoint3$surv[2] + 0.03, label = paste0(signif(timepoint3$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 12, y = timepoint4$surv[1] - 0.08, label = paste0(signif(timepoint4$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 12, y = timepoint4$surv[2] + 0.03, label = paste0(signif(timepoint4$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 15, y = timepoint5$surv[1] - 0.08, label = paste0(signif(timepoint5$surv[1]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 15, y = timepoint5$surv[2] + 0.03, label = paste0(signif(timepoint5$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = -0.7, y = 0,
           label = paste0("Elra vs tec\n",
                          " Median: ", signif(med_time$median[1], 3), " vs ", signif(med_time$median[2], 3), " months\n",
                          " HR: ", sprintf("%.2f", HR), " (95% CI ", sprintf("%.2f", CI_lower), "-", sprintf("%.2f", CI_upper),", p=",pval, ")"),
           color = "black", hjust = 0, vjust = 0, size = 3)

DOR

6.3 ***FIGURE 3: PFS/OS/DOR, stratified by trial elibility

FIGURE3 <- wrap_plots(
  (OS$plot / OS$table) + plot_layout(heights = c(4, 0.9)),
  (PFS$plot / PFS$table) + plot_layout(heights = c(4, 0.9)),
  (DOR$plot / DOR$table) + plot_layout(heights = c(4, 0.9)),
  ncol = 1  # Force vertical stacking
)

ggsave("svg/FIGURE3.svg", FIGURE3, device = "svg",
       width = 6.5,
       height = 8.5,
       units = c("in"))

ggsave("jpeg/FIGURE3.jpeg", FIGURE3, device = "jpeg",
       width = 6.5,
       height = 8.5,
       dpi = 600,
       units = c("in"))

ggsave("tiff/FIGURE3.tiff", FIGURE3, device = "tiff",
       width = 6.5,
       height = 8.5,
       dpi = 600,
       units = c("in"),
       compression = "lzw")

knitr::include_graphics("jpeg/FIGURE3.jpeg")

6.3.1 Schoenfeld residuals and PH test

data1 <- IPTW.data %>%
  transmute(
    site,
    patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate,
    elra_ind = ifelse(bs_ab_name == "elra", 1, 0),  # 1 = elra, 0 = tec
    Relapse.or.death = ifelse(progression_yes_no | death_yes_no,1,0),
    Death = ifelse(death_yes_no,1,0),
    date_of_bs_ab_step_up_d1,
    date_of_pd,
    date_of_last_contact,
    date_progression_or_death = coalesce(date_of_pd, date_of_last_contact),
    Days.to.death.or.DLC = as.numeric(date_of_last_contact - date_of_bs_ab_step_up_d1),
    Days.to.death.or.relapse.or.DLC = as.numeric(date_progression_or_death  - date_of_bs_ab_step_up_d1),
    trial_eligible = ifelse(eligible_for_corresponding_bs_ab_trial_on_first_dose_use_trial_elig_sheet_to_help_you ==1,1,0),
    weights
  )

## ---- Cox models ----
cox_OS <- coxph(Surv(Days.to.death.or.DLC/30, Death) ~ elra_ind, weights = weights, data = data1)

cox_PFS <- coxph(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ elra_ind, weights = weights, data = data1)

## ---- PH diagnostics ----
zph_OS  <- cox.zph(cox_OS)
zph_PFS <- cox.zph(cox_PFS)

# Global PH p-values (same as elra_ind right now)
p_OS  <- signif(zph_OS$table["GLOBAL",  "p"], 3)
p_PFS <- signif(zph_PFS$table["GLOBAL", "p"], 3)

par(mfrow = c(1, 2))

plot(
  zph_OS,
  main = paste0(
    "Schoenfeld residuals – OS\n",
    "Global PH test: p = ", p_OS
  )
)

plot(
  zph_PFS,
  main = paste0(
    "Schoenfeld residuals – PFS\n",
    "Global PH test: p = ", p_PFS
  )
)

par(mfrow = c(1, 1))

6.4 High risk subsets

6.4.1 True EMD

data1 <- IPTW.data %>%
  filter(EMD == 1) %>%
  transmute(
    site,
    patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate,
    bs_ab_name,
    Relapse.or.death = ifelse(progression_yes_no | death_yes_no,1,0),
    Death = ifelse(death_yes_no,1,0),
    date_of_bs_ab_step_up_d1,
    date_of_pd,
    date_of_last_contact,
    date_progression_or_death = coalesce(date_of_pd, date_of_last_contact),
    Days.to.death.or.DLC = as.numeric(date_of_last_contact - date_of_bs_ab_step_up_d1),
    Days.to.death.or.relapse.or.DLC = as.numeric(date_progression_or_death  - date_of_bs_ab_step_up_d1),
    trial_eligible = ifelse(eligible_for_corresponding_bs_ab_trial_on_first_dose_use_trial_elig_sheet_to_help_you ==1,1,0),
    weights
  )

data_dor <- IPTW.data %>%
  filter(best_response %in% c("sCR","CR","VGPR","PR")) %>%
  filter(EMD == 1) %>%
  transmute(
    weights,
    bs_ab_name,
    day_30_response,
    day_90_response,
    best_response,
    Relapse.or.death = ifelse(progression_yes_no | death_yes_no,1,0),
    Death = ifelse(death_yes_no,1,0),
    date_of_pd,
    date_of_last_contact,
    date_of_bs_ab_step_up_d1,
    date_progression_or_death = coalesce(date_of_pd, date_of_last_contact),
    day_30_response,
    day_90_response,
    best_response,
    time_to_best_response = case_when(
      best_response == day_30_response ~ 30,
      best_response == day_90_response ~ 90,
      .default = NA
    ),
    trial_eligible = ifelse(eligible_for_corresponding_bs_ab_trial_on_first_dose_use_trial_elig_sheet_to_help_you ==1,1,0)
  ) %>%
  filter(!is.na(time_to_best_response)) %>%
  mutate(
    Days.to.death.or.DLC = as.numeric(date_of_last_contact - date_of_bs_ab_step_up_d1) - time_to_best_response,
    Days.to.death.or.relapse.or.DLC = as.numeric(date_progression_or_death  - date_of_bs_ab_step_up_d1)- time_to_best_response,
  ) %>%
  filter(Days.to.death.or.DLC >0 & Days.to.death.or.relapse.or.DLC >0)

km_OS <- survfit(Surv(Days.to.death.or.DLC/30, Death) ~ bs_ab_name, weights = weights, data = data1)
km_OS
## Call: survfit(formula = Surv(Days.to.death.or.DLC/30, Death) ~ bs_ab_name, 
##     data = data1, weights = weights)
## 
##    4 observations deleted due to missingness 
##                 records    n events median 0.95LCL 0.95UCL
## bs_ab_name=elra      24 18.2   4.86     NA    5.23      NA
## bs_ab_name=tec       32 37.4  21.05   14.1   10.43      NA
med_time <- surv_median(km_OS)

cox_model <- coxph(Surv(Days.to.death.or.DLC/30, Death) ~ ifelse(bs_ab_name == "elra",1,0), weights = weights, data = data1)
cox_summary <- summary(cox_model)

HR <- round(cox_summary$coefficients[1, "exp(coef)"], 2)
CI_lower <- round(cox_summary$conf.int[1, "lower .95"], 2)
CI_upper <- round(cox_summary$conf.int[1, "upper .95"], 2)
pval <- signif(cox_summary$coefficients[1, "Pr(>|z|)"], 2)

## 1-year timepoint
timepoint1 <- summary(km_OS,times=c(3))
timepoint2 <- summary(km_OS,times=c(6))
timepoint3 <- summary(km_OS,times=c(9))
timepoint4 <- summary(km_OS,times=c(12))
timepoint5 <- summary(km_OS,times=c(15))

OS <- ggsurvplot(km_OS,
           pval=FALSE,
           conf.int = FALSE,
           risk.table = TRUE, # Add risk table
           fontsize = 4,
           risk.table.col = "strata", # Change risk table color by groups
           tables.height = 0.3,
           tables.theme = theme_cleantable() +
             theme(
               axis.text.y = element_text(size = 9),
               plot.title = element_text(size = 10)
               ),
           xlab="Time (months)", ylab = "Overall survival",
           xlim = c(0,18), break.x.by = c(3),
           ylim = c(0,1), break.y.by = c(0.25),
           legend = "none",
           surv.scale="percent",
           font.main = c(10, "plain", "black"),
           font.x = c(10, "plain", "black"),
           font.y = c(10, "plain", "black"),
           font.caption = c(10, "plain", "black"),
           font.tickslab = c(10, "plain", "black"),
           legend.labs = c("Elra", "Tec"),
           title = "A)",
           size = 0.5
) 

OS$plot <- OS$plot + 
  geom_vline(xintercept = c(3, 6, 9, 12, 15), linetype = "dotted", color = "gray50") +
  annotate("text", x = 3, y = timepoint1$surv[1] - 0.1, label = paste0(signif(timepoint1$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 3, y = timepoint1$surv[2] + 0.03, label = paste0(signif(timepoint1$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 6, y = timepoint2$surv[1] - 0.08, label = paste0(signif(timepoint2$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 6, y = timepoint2$surv[2] + 0.03, label = paste0(signif(timepoint2$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 9, y = timepoint3$surv[1] - 0.08, label = paste0(signif(timepoint3$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 9, y = timepoint3$surv[2] + 0.03, label = paste0(signif(timepoint3$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 12, y = timepoint4$surv[1] - 0.08, label = paste0(signif(timepoint4$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 12, y = timepoint4$surv[2] + 0.03, label = paste0(signif(timepoint4$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 15, y = timepoint5$surv[1] - 0.08, label = paste0(signif(timepoint5$surv[1]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 15, y = timepoint5$surv[2] + 0.03, label = paste0(signif(timepoint5$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = -0.7, y = 0,
           label = paste0("Elra vs tec\n",
                          " Median: ", signif(med_time$median[1], 3), " vs ", signif(med_time$median[2], 3), " months\n",
                          " HR: ", sprintf("%.2f", HR), " (95% CI ", sprintf("%.2f", CI_lower), "-", sprintf("%.2f", CI_upper),", p=",pval, ")"),
           color = "black", hjust = 0, vjust = 0, size = 3)

OS

km_PFS <- survfit(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ bs_ab_name, weights = weights, data = data1)
km_PFS
## Call: survfit(formula = Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ 
##     bs_ab_name, data = data1, weights = weights)
## 
##    1 observation deleted due to missingness 
##                 records    n events median 0.95LCL 0.95UCL
## bs_ab_name=elra      25 18.8   8.27   6.27    2.57      NA
## bs_ab_name=tec       34 39.5  27.45   7.93    4.67    12.5
med_time <- surv_median(km_PFS)

cox_model <- coxph(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ ifelse(bs_ab_name == "elra",1,0), weights = weights, data = data1)
cox_summary <- summary(cox_model)

HR <- round(cox_summary$coefficients[1, "exp(coef)"], 2)
CI_lower <- round(cox_summary$conf.int[1, "lower .95"], 2)
CI_upper <- round(cox_summary$conf.int[1, "upper .95"], 2)
pval <- signif(cox_summary$coefficients[1, "Pr(>|z|)"], 2)

## Timepoints
timepoint1 <- summary(km_PFS,times=c(3))
timepoint2 <- summary(km_PFS,times=c(6))
timepoint3 <- summary(km_PFS,times=c(9))
timepoint4 <- summary(km_PFS,times=c(12))
timepoint5 <- summary(km_PFS,times=c(15))


PFS <- ggsurvplot(km_PFS,
           pval=FALSE,
           conf.int = FALSE,
           risk.table = TRUE, # Add risk table
           fontsize = 4,
           risk.table.col = "strata", # Change risk table color by groups
           tables.height = 0.3,
           tables.theme = theme_cleantable() +
             theme(
               axis.text.y = element_text(size = 9),
               plot.title = element_text(size = 10)
               ),
           xlab="Time (months)", ylab = "Progression-free survival",
           xlim = c(0,18), break.x.by = c(3),
           ylim = c(0,1), break.y.by = c(0.25),
           legend = "none",
           surv.scale="percent",
           font.main = c(10, "plain", "black"),
           font.x = c(10, "plain", "black"),
           font.y = c(10, "plain", "black"),
           font.caption = c(10, "plain", "black"),
           font.tickslab = c(10, "plain", "black"),
           legend.labs = c("Elra", "Tec"),
           title = "B)",
           size = 0.5
)

PFS$plot <- PFS$plot + 
  geom_vline(xintercept = c(3, 6, 9, 12, 15), linetype = "dotted", color = "gray50") +
  annotate("text", x = 3, y = timepoint1$surv[1] - 0.08, label = paste0(signif(timepoint1$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 3, y = timepoint1$surv[2] + 0.03, label = paste0(signif(timepoint1$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 6, y = timepoint2$surv[1] - 0.08, label = paste0(signif(timepoint2$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 6, y = timepoint2$surv[2] + 0.03, label = paste0(signif(timepoint2$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 9, y = timepoint3$surv[1] - 0.08, label = paste0(signif(timepoint3$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 9, y = timepoint3$surv[2] + 0.03, label = paste0(signif(timepoint3$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 12, y = timepoint4$surv[1] - 0.08, label = paste0(signif(timepoint4$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 12, y = timepoint4$surv[2] + 0.03, label = paste0(signif(timepoint4$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 15, y = timepoint5$surv[1] - 0.08, label = paste0(signif(timepoint5$surv[1]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 15, y = timepoint5$surv[2] + 0.03, label = paste0(signif(timepoint5$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = -0.7, y = 0,
           label = paste0("Elra vs tec\n",
                          " Median: ", signif(med_time$median[1], 3), " vs ", signif(med_time$median[2], 3), " months\n",
                          " HR: ", sprintf("%.2f", HR), " (95% CI ", sprintf("%.2f", CI_lower), "-", sprintf("%.2f", CI_upper),", p=",pval, ")"),
           color = "black", hjust = 0, vjust = 0, size = 3)

PFS

km_DOR <- survfit(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ bs_ab_name, weights = weights, data = data_dor)
km_DOR
## Call: survfit(formula = Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ 
##     bs_ab_name, data = data_dor, weights = weights)
## 
##    1 observation deleted due to missingness 
##                 records     n events median 0.95LCL 0.95UCL
## bs_ab_name=elra      11  8.76   2.39     NA    5.27      NA
## bs_ab_name=tec       15 18.05   8.69   6.93    5.63      NA
med_time <- surv_median(km_DOR)

cox_model <- coxph(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ ifelse(bs_ab_name == "elra",1,0), weights = weights, data = data_dor)
cox_summary <- summary(cox_model)

HR <- round(cox_summary$coefficients[1, "exp(coef)"], 2)
CI_lower <- round(cox_summary$conf.int[1, "lower .95"], 2)
CI_upper <- round(cox_summary$conf.int[1, "upper .95"], 2)
pval <- signif(cox_summary$coefficients[1, "Pr(>|z|)"], 2)

## Specific timepoints
timepoint1 <- summary(km_DOR,times=c(3))
timepoint2 <- summary(km_DOR,times=c(6))
timepoint3 <- summary(km_DOR,times=c(9))
timepoint4 <- summary(km_DOR,times=c(12))
timepoint5 <- summary(km_DOR,times=c(15))


DOR <- ggsurvplot(km_DOR,
           pval=FALSE,
           conf.int = FALSE,
           risk.table = TRUE, # Add risk table
           fontsize = 4,
           risk.table.col = "strata", # Change risk table color by groups
           tables.height = 0.3,
           #risk.table.fontsize = 4,
           tables.theme = theme_cleantable() +
             theme(
               axis.text.y = element_text(size = 9),
               plot.title = element_text(size = 10)
               ),
           xlab="Time (months)", ylab = "Duration of response",
           xlim = c(0,18), break.x.by = c(3),
           ylim = c(0,1), break.y.by = c(0.25),
           legend = "none",
           surv.scale="percent",
           font.main = c(10, "plain", "black"),
           font.x = c(10, "plain", "black"),
           font.y = c(10, "plain", "black"),
           font.caption = c(10, "plain", "black"),
           font.tickslab = c(10, "plain", "black"),
           legend.labs = c("Elra", "Tec"),
           title = "C)",
           size = 0.5
)

DOR$plot <- DOR$plot + 
  geom_vline(xintercept = c(3, 6, 9, 12, 15), linetype = "dotted", color = "gray50") +
  annotate("text", x = 3, y = timepoint1$surv[1] -0.08, label = paste0(signif(timepoint1$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 3, y = timepoint1$surv[2] + 0.03, label = paste0(signif(timepoint1$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 6, y = timepoint2$surv[1] - 0.08, label = paste0(signif(timepoint2$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 6, y = timepoint2$surv[2] + 0.03, label = paste0(signif(timepoint2$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 9, y = timepoint3$surv[1] - 0.08, label = paste0(signif(timepoint3$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 9, y = timepoint3$surv[2] + 0.03, label = paste0(signif(timepoint3$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 12, y = timepoint4$surv[1] - 0.08, label = paste0(signif(timepoint4$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 12, y = timepoint4$surv[2] + 0.03, label = paste0(signif(timepoint4$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 15, y = timepoint5$surv[1] - 0.08, label = paste0(signif(timepoint5$surv[1]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 15, y = timepoint5$surv[2] + 0.03, label = paste0(signif(timepoint5$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = -0.7, y = 0,
           label = paste0("Elra vs tec\n",
                          " Median: ", signif(med_time$median[1], 3), " vs ", signif(med_time$median[2], 3), " months\n",
                          " HR: ", sprintf("%.2f", HR), " (95% CI ", sprintf("%.2f", CI_lower), "-", sprintf("%.2f", CI_upper),", p=",pval, ")"),
           color = "black", hjust = 0, vjust = 0, size = 3)

DOR

6.4.2 ECOG 2+

data1 <- IPTW.data %>%
  filter(ecog_ps_at_bs_ab_initiation_0_5 >= 2) %>%
  transmute(
    site,
    patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate,
    bs_ab_name,
    Relapse.or.death = ifelse(progression_yes_no | death_yes_no,1,0),
    Death = ifelse(death_yes_no,1,0),
    date_of_bs_ab_step_up_d1,
    date_of_pd,
    date_of_last_contact,
    date_progression_or_death = coalesce(date_of_pd, date_of_last_contact),
    Days.to.death.or.DLC = as.numeric(date_of_last_contact - date_of_bs_ab_step_up_d1),
    Days.to.death.or.relapse.or.DLC = as.numeric(date_progression_or_death  - date_of_bs_ab_step_up_d1),
    trial_eligible = ifelse(eligible_for_corresponding_bs_ab_trial_on_first_dose_use_trial_elig_sheet_to_help_you ==1,1,0),
    weights
  )

data_dor <- IPTW.data %>%
  filter(best_response %in% c("sCR","CR","VGPR","PR")) %>%
  filter(ecog_ps_at_bs_ab_initiation_0_5 >= 2) %>%
  transmute(
    weights,
    bs_ab_name,
    day_30_response,
    day_90_response,
    best_response,
    Relapse.or.death = ifelse(progression_yes_no | death_yes_no,1,0),
    Death = ifelse(death_yes_no,1,0),
    date_of_pd,
    date_of_last_contact,
    date_of_bs_ab_step_up_d1,
    date_progression_or_death = coalesce(date_of_pd, date_of_last_contact),
    day_30_response,
    day_90_response,
    best_response,
    time_to_best_response = case_when(
      best_response == day_30_response ~ 30,
      best_response == day_90_response ~ 90,
      .default = NA
    ),
    trial_eligible = ifelse(eligible_for_corresponding_bs_ab_trial_on_first_dose_use_trial_elig_sheet_to_help_you ==1,1,0)
  ) %>%
  filter(!is.na(time_to_best_response)) %>%
  mutate(
    Days.to.death.or.DLC = as.numeric(date_of_last_contact - date_of_bs_ab_step_up_d1) - time_to_best_response,
    Days.to.death.or.relapse.or.DLC = as.numeric(date_progression_or_death  - date_of_bs_ab_step_up_d1)- time_to_best_response,
  ) %>%
  filter(Days.to.death.or.DLC >0 & Days.to.death.or.relapse.or.DLC >0)

km_OS <- survfit(Surv(Days.to.death.or.DLC/30, Death) ~ bs_ab_name, weights = weights, data = data1)
km_OS
## Call: survfit(formula = Surv(Days.to.death.or.DLC/30, Death) ~ bs_ab_name, 
##     data = data1, weights = weights)
## 
##    4 observations deleted due to missingness 
##                 records    n events median 0.95LCL 0.95UCL
## bs_ab_name=elra      42 37.9   18.8   4.87    2.27      NA
## bs_ab_name=tec       74 75.8   45.8   9.57    5.27    18.8
med_time <- surv_median(km_OS)

cox_model <- coxph(Surv(Days.to.death.or.DLC/30, Death) ~ ifelse(bs_ab_name == "elra",1,0), weights = weights, data = data1)
cox_summary <- summary(cox_model)

HR <- round(cox_summary$coefficients[1, "exp(coef)"], 2)
CI_lower <- round(cox_summary$conf.int[1, "lower .95"], 2)
CI_upper <- round(cox_summary$conf.int[1, "upper .95"], 2)
pval <- signif(cox_summary$coefficients[1, "Pr(>|z|)"], 2)

## 1-year timepoint
timepoint1 <- summary(km_OS,times=c(3))
timepoint2 <- summary(km_OS,times=c(6))
timepoint3 <- summary(km_OS,times=c(9))
timepoint4 <- summary(km_OS,times=c(12))
timepoint5 <- summary(km_OS,times=c(15))

OS <- ggsurvplot(km_OS,
           pval=FALSE,
           conf.int = FALSE,
           risk.table = TRUE, # Add risk table
           fontsize = 4,
           risk.table.col = "strata", # Change risk table color by groups
           tables.height = 0.3,
           tables.theme = theme_cleantable() +
             theme(
               axis.text.y = element_text(size = 9),
               plot.title = element_text(size = 10)
               ),
           xlab="Time (months)", ylab = "Overall survival",
           xlim = c(0,18), break.x.by = c(3),
           ylim = c(0,1), break.y.by = c(0.25),
           legend = "none",
           surv.scale="percent",
           font.main = c(10, "plain", "black"),
           font.x = c(10, "plain", "black"),
           font.y = c(10, "plain", "black"),
           font.caption = c(10, "plain", "black"),
           font.tickslab = c(10, "plain", "black"),
           legend.labs = c("Elra", "Tec"),
           title = "A)",
           size = 0.5
) 

OS$plot <- OS$plot + 
  geom_vline(xintercept = c(3, 6, 9, 12, 15), linetype = "dotted", color = "gray50") +
  annotate("text", x = 3, y = timepoint1$surv[1] - 0.1, label = paste0(signif(timepoint1$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 3, y = timepoint1$surv[2] + 0.03, label = paste0(signif(timepoint1$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 6, y = timepoint2$surv[1] - 0.08, label = paste0(signif(timepoint2$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 6, y = timepoint2$surv[2] + 0.03, label = paste0(signif(timepoint2$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 9, y = timepoint3$surv[1] - 0.08, label = paste0(signif(timepoint3$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 9, y = timepoint3$surv[2] + 0.03, label = paste0(signif(timepoint3$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 12, y = timepoint4$surv[1] - 0.08, label = paste0(signif(timepoint4$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 12, y = timepoint4$surv[2] + 0.03, label = paste0(signif(timepoint4$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 15, y = timepoint5$surv[1] - 0.08, label = paste0(signif(timepoint5$surv[1]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 15, y = timepoint5$surv[2] + 0.03, label = paste0(signif(timepoint5$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = -0.7, y = 0,
           label = paste0("Elra vs tec\n",
                          " Median: ", signif(med_time$median[1], 3), " vs ", signif(med_time$median[2], 3), " months\n",
                          " HR: ", sprintf("%.2f", HR), " (95% CI ", sprintf("%.2f", CI_lower), "-", sprintf("%.2f", CI_upper),", p=",pval, ")"),
           color = "black", hjust = 0, vjust = 0, size = 3)

OS

km_PFS <- survfit(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ bs_ab_name, weights = weights, data = data1)
km_PFS
## Call: survfit(formula = Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ 
##     bs_ab_name, data = data1, weights = weights)
## 
##    2 observations deleted due to missingness 
##                 records    n events median 0.95LCL 0.95UCL
## bs_ab_name=elra      43 38.6   27.9   1.63    1.07   10.27
## bs_ab_name=tec       75 76.7   56.7   4.43    2.10    8.73
med_time <- surv_median(km_PFS)

cox_model <- coxph(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ ifelse(bs_ab_name == "elra",1,0), weights = weights, data = data1)
cox_summary <- summary(cox_model)

HR <- round(cox_summary$coefficients[1, "exp(coef)"], 2)
CI_lower <- round(cox_summary$conf.int[1, "lower .95"], 2)
CI_upper <- round(cox_summary$conf.int[1, "upper .95"], 2)
pval <- signif(cox_summary$coefficients[1, "Pr(>|z|)"], 2)

## Timepoints
timepoint1 <- summary(km_PFS,times=c(3))
timepoint2 <- summary(km_PFS,times=c(6))
timepoint3 <- summary(km_PFS,times=c(9))
timepoint4 <- summary(km_PFS,times=c(12))
timepoint5 <- summary(km_PFS,times=c(15))


PFS <- ggsurvplot(km_PFS,
           pval=FALSE,
           conf.int = FALSE,
           risk.table = TRUE, # Add risk table
           fontsize = 4,
           risk.table.col = "strata", # Change risk table color by groups
           tables.height = 0.3,
           tables.theme = theme_cleantable() +
             theme(
               axis.text.y = element_text(size = 9),
               plot.title = element_text(size = 10)
               ),
           xlab="Time (months)", ylab = "Progression-free survival",
           xlim = c(0,18), break.x.by = c(3),
           ylim = c(0,1), break.y.by = c(0.25),
           legend = "none",
           surv.scale="percent",
           font.main = c(10, "plain", "black"),
           font.x = c(10, "plain", "black"),
           font.y = c(10, "plain", "black"),
           font.caption = c(10, "plain", "black"),
           font.tickslab = c(10, "plain", "black"),
           legend.labs = c("Elra", "Tec"),
           title = "B)",
           size = 0.5
)

PFS$plot <- PFS$plot + 
  geom_vline(xintercept = c(3, 6, 9, 12, 15), linetype = "dotted", color = "gray50") +
  annotate("text", x = 3, y = timepoint1$surv[1] - 0.08, label = paste0(signif(timepoint1$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 3, y = timepoint1$surv[2] + 0.03, label = paste0(signif(timepoint1$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 6, y = timepoint2$surv[1] - 0.08, label = paste0(signif(timepoint2$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 6, y = timepoint2$surv[2] + 0.03, label = paste0(signif(timepoint2$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 9, y = timepoint3$surv[1] - 0.08, label = paste0(signif(timepoint3$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 9, y = timepoint3$surv[2] + 0.03, label = paste0(signif(timepoint3$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 12, y = timepoint4$surv[1] - 0.08, label = paste0(signif(timepoint4$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 12, y = timepoint4$surv[2] + 0.03, label = paste0(signif(timepoint4$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 15, y = timepoint5$surv[1] - 0.08, label = paste0(signif(timepoint5$surv[1]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 15, y = timepoint5$surv[2] + 0.03, label = paste0(signif(timepoint5$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = -0.7, y = 0,
           label = paste0("Elra vs tec\n",
                          " Median: ", signif(med_time$median[1], 3), " vs ", signif(med_time$median[2], 3), " months\n",
                          " HR: ", sprintf("%.2f", HR), " (95% CI ", sprintf("%.2f", CI_lower), "-", sprintf("%.2f", CI_upper),", p=",pval, ")"),
           color = "black", hjust = 0, vjust = 0, size = 3)

PFS

km_DOR <- survfit(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ bs_ab_name, weights = weights, data = data_dor)
km_DOR
## Call: survfit(formula = Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ 
##     bs_ab_name, data = data_dor, weights = weights)
## 
##                 records    n events median 0.95LCL 0.95UCL
## bs_ab_name=elra      15 12.9   6.02   11.7     1.0      NA
## bs_ab_name=tec       30 31.3  16.87   13.3     8.3      NA
med_time <- surv_median(km_DOR)

cox_model <- coxph(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ ifelse(bs_ab_name == "elra",1,0), weights = weights, data = data_dor)
cox_summary <- summary(cox_model)

HR <- round(cox_summary$coefficients[1, "exp(coef)"], 2)
CI_lower <- round(cox_summary$conf.int[1, "lower .95"], 2)
CI_upper <- round(cox_summary$conf.int[1, "upper .95"], 2)
pval <- signif(cox_summary$coefficients[1, "Pr(>|z|)"], 2)

## Specific timepoints
timepoint1 <- summary(km_DOR,times=c(3))
timepoint2 <- summary(km_DOR,times=c(6))
timepoint3 <- summary(km_DOR,times=c(9))
timepoint4 <- summary(km_DOR,times=c(12))
timepoint5 <- summary(km_DOR,times=c(15))


DOR <- ggsurvplot(km_DOR,
           pval=FALSE,
           conf.int = FALSE,
           risk.table = TRUE, # Add risk table
           fontsize = 4,
           risk.table.col = "strata", # Change risk table color by groups
           tables.height = 0.3,
           #risk.table.fontsize = 4,
           tables.theme = theme_cleantable() +
             theme(
               axis.text.y = element_text(size = 9),
               plot.title = element_text(size = 10)
               ),
           xlab="Time (months)", ylab = "Duration of response",
           xlim = c(0,18), break.x.by = c(3),
           ylim = c(0,1), break.y.by = c(0.25),
           legend = "none",
           surv.scale="percent",
           font.main = c(10, "plain", "black"),
           font.x = c(10, "plain", "black"),
           font.y = c(10, "plain", "black"),
           font.caption = c(10, "plain", "black"),
           font.tickslab = c(10, "plain", "black"),
           legend.labs = c("Elra", "Tec"),
           title = "C)",
           size = 0.5
)

DOR$plot <- DOR$plot + 
  geom_vline(xintercept = c(3, 6, 9, 12, 15), linetype = "dotted", color = "gray50") +
  annotate("text", x = 3, y = timepoint1$surv[1] -0.08, label = paste0(signif(timepoint1$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 3, y = timepoint1$surv[2] + 0.03, label = paste0(signif(timepoint1$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 6, y = timepoint2$surv[1] - 0.08, label = paste0(signif(timepoint2$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 6, y = timepoint2$surv[2] + 0.03, label = paste0(signif(timepoint2$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 9, y = timepoint3$surv[1] - 0.08, label = paste0(signif(timepoint3$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 9, y = timepoint3$surv[2] + 0.03, label = paste0(signif(timepoint3$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 12, y = timepoint4$surv[1] - 0.08, label = paste0(signif(timepoint4$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 12, y = timepoint4$surv[2] + 0.03, label = paste0(signif(timepoint4$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 15, y = timepoint5$surv[1] - 0.08, label = paste0(signif(timepoint5$surv[1]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 15, y = timepoint5$surv[2] + 0.03, label = paste0(signif(timepoint5$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = -0.7, y = 0,
           label = paste0("Elra vs tec\n",
                          " Median: ", signif(med_time$median[1], 3), " vs ", signif(med_time$median[2], 3), " months\n",
                          " HR: ", sprintf("%.2f", HR), " (95% CI ", sprintf("%.2f", CI_lower), "-", sprintf("%.2f", CI_upper),", p=",pval, ")"),
           color = "black", hjust = 0, vjust = 0, size = 3)

DOR

6.4.3 Hgb <10 g/dL

data1 <- IPTW.data %>%
  filter(baseline_hgb <10) %>%
  transmute(
    site,
    patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate,
    bs_ab_name,
    Relapse.or.death = ifelse(progression_yes_no | death_yes_no,1,0),
    Death = ifelse(death_yes_no,1,0),
    date_of_bs_ab_step_up_d1,
    date_of_pd,
    date_of_last_contact,
    date_progression_or_death = coalesce(date_of_pd, date_of_last_contact),
    Days.to.death.or.DLC = as.numeric(date_of_last_contact - date_of_bs_ab_step_up_d1),
    Days.to.death.or.relapse.or.DLC = as.numeric(date_progression_or_death  - date_of_bs_ab_step_up_d1),
    trial_eligible = ifelse(eligible_for_corresponding_bs_ab_trial_on_first_dose_use_trial_elig_sheet_to_help_you ==1,1,0),
    weights
  )

data_dor <- IPTW.data %>%
  filter(best_response %in% c("sCR","CR","VGPR","PR")) %>%
  filter(baseline_hgb <10) %>%
  transmute(
    weights,
    bs_ab_name,
    day_30_response,
    day_90_response,
    best_response,
    Relapse.or.death = ifelse(progression_yes_no | death_yes_no,1,0),
    Death = ifelse(death_yes_no,1,0),
    date_of_pd,
    date_of_last_contact,
    date_of_bs_ab_step_up_d1,
    date_progression_or_death = coalesce(date_of_pd, date_of_last_contact),
    day_30_response,
    day_90_response,
    best_response,
    time_to_best_response = case_when(
      best_response == day_30_response ~ 30,
      best_response == day_90_response ~ 90,
      .default = NA
    ),
    trial_eligible = ifelse(eligible_for_corresponding_bs_ab_trial_on_first_dose_use_trial_elig_sheet_to_help_you ==1,1,0)
  ) %>%
  filter(!is.na(time_to_best_response)) %>%
  mutate(
    Days.to.death.or.DLC = as.numeric(date_of_last_contact - date_of_bs_ab_step_up_d1) - time_to_best_response,
    Days.to.death.or.relapse.or.DLC = as.numeric(date_progression_or_death  - date_of_bs_ab_step_up_d1)- time_to_best_response,
  ) %>%
  filter(Days.to.death.or.DLC >0 & Days.to.death.or.relapse.or.DLC >0)

km_OS <- survfit(Surv(Days.to.death.or.DLC/30, Death) ~ bs_ab_name, weights = weights, data = data1)
km_OS
## Call: survfit(formula = Surv(Days.to.death.or.DLC/30, Death) ~ bs_ab_name, 
##     data = data1, weights = weights)
## 
##    8 observations deleted due to missingness 
##                 records     n events median 0.95LCL 0.95UCL
## bs_ab_name=elra      55  63.3   31.7   4.87    3.17      NA
## bs_ab_name=tec      134 125.3   71.5  12.83    9.57    18.6
med_time <- surv_median(km_OS)

cox_model <- coxph(Surv(Days.to.death.or.DLC/30, Death) ~ ifelse(bs_ab_name == "elra",1,0), weights = weights, data = data1)
cox_summary <- summary(cox_model)

HR <- round(cox_summary$coefficients[1, "exp(coef)"], 2)
CI_lower <- round(cox_summary$conf.int[1, "lower .95"], 2)
CI_upper <- round(cox_summary$conf.int[1, "upper .95"], 2)
pval <- signif(cox_summary$coefficients[1, "Pr(>|z|)"], 2)

## 1-year timepoint
timepoint1 <- summary(km_OS,times=c(3))
timepoint2 <- summary(km_OS,times=c(6))
timepoint3 <- summary(km_OS,times=c(9))
timepoint4 <- summary(km_OS,times=c(12))
timepoint5 <- summary(km_OS,times=c(15))

OS <- ggsurvplot(km_OS,
           pval=FALSE,
           conf.int = FALSE,
           risk.table = TRUE, # Add risk table
           fontsize = 4,
           risk.table.col = "strata", # Change risk table color by groups
           tables.height = 0.3,
           tables.theme = theme_cleantable() +
             theme(
               axis.text.y = element_text(size = 9),
               plot.title = element_text(size = 10)
               ),
           xlab="Time (months)", ylab = "Overall survival",
           xlim = c(0,18), break.x.by = c(3),
           ylim = c(0,1), break.y.by = c(0.25),
           legend = "none",
           surv.scale="percent",
           font.main = c(10, "plain", "black"),
           font.x = c(10, "plain", "black"),
           font.y = c(10, "plain", "black"),
           font.caption = c(10, "plain", "black"),
           font.tickslab = c(10, "plain", "black"),
           legend.labs = c("Elra", "Tec"),
           title = "A)",
           size = 0.5
) 

OS$plot <- OS$plot + 
  geom_vline(xintercept = c(3, 6, 9, 12, 15), linetype = "dotted", color = "gray50") +
  annotate("text", x = 3, y = timepoint1$surv[1] - 0.1, label = paste0(signif(timepoint1$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 3, y = timepoint1$surv[2] + 0.03, label = paste0(signif(timepoint1$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 6, y = timepoint2$surv[1] - 0.08, label = paste0(signif(timepoint2$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 6, y = timepoint2$surv[2] + 0.03, label = paste0(signif(timepoint2$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 9, y = timepoint3$surv[1] - 0.08, label = paste0(signif(timepoint3$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 9, y = timepoint3$surv[2] + 0.03, label = paste0(signif(timepoint3$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 12, y = timepoint4$surv[1] - 0.08, label = paste0(signif(timepoint4$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 12, y = timepoint4$surv[2] + 0.03, label = paste0(signif(timepoint4$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 15, y = timepoint5$surv[1] - 0.08, label = paste0(signif(timepoint5$surv[1]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 15, y = timepoint5$surv[2] + 0.03, label = paste0(signif(timepoint5$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = -0.7, y = 0,
           label = paste0("Elra vs tec\n",
                          " Median: ", signif(med_time$median[1], 3), " vs ", signif(med_time$median[2], 3), " months\n",
                          " HR: ", sprintf("%.2f", HR), " (95% CI ", sprintf("%.2f", CI_lower), "-", sprintf("%.2f", CI_upper),", p=",pval, ")"),
           color = "black", hjust = 0, vjust = 0, size = 3)

OS

km_PFS <- survfit(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ bs_ab_name, weights = weights, data = data1)
km_PFS
## Call: survfit(formula = Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ 
##     bs_ab_name, data = data1, weights = weights)
## 
##    4 observations deleted due to missingness 
##                 records     n events median 0.95LCL 0.95UCL
## bs_ab_name=elra      56  63.9   43.7   1.63    1.10    6.27
## bs_ab_name=tec      137 128.3   97.8   3.70    2.17    7.40
med_time <- surv_median(km_PFS)

cox_model <- coxph(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ ifelse(bs_ab_name == "elra",1,0), weights = weights, data = data1)
cox_summary <- summary(cox_model)

HR <- round(cox_summary$coefficients[1, "exp(coef)"], 2)
CI_lower <- round(cox_summary$conf.int[1, "lower .95"], 2)
CI_upper <- round(cox_summary$conf.int[1, "upper .95"], 2)
pval <- signif(cox_summary$coefficients[1, "Pr(>|z|)"], 2)

## Timepoints
timepoint1 <- summary(km_PFS,times=c(3))
timepoint2 <- summary(km_PFS,times=c(6))
timepoint3 <- summary(km_PFS,times=c(9))
timepoint4 <- summary(km_PFS,times=c(12))
timepoint5 <- summary(km_PFS,times=c(15))


PFS <- ggsurvplot(km_PFS,
           pval=FALSE,
           conf.int = FALSE,
           risk.table = TRUE, # Add risk table
           fontsize = 4,
           risk.table.col = "strata", # Change risk table color by groups
           tables.height = 0.3,
           tables.theme = theme_cleantable() +
             theme(
               axis.text.y = element_text(size = 9),
               plot.title = element_text(size = 10)
               ),
           xlab="Time (months)", ylab = "Progression-free survival",
           xlim = c(0,18), break.x.by = c(3),
           ylim = c(0,1), break.y.by = c(0.25),
           legend = "none",
           surv.scale="percent",
           font.main = c(10, "plain", "black"),
           font.x = c(10, "plain", "black"),
           font.y = c(10, "plain", "black"),
           font.caption = c(10, "plain", "black"),
           font.tickslab = c(10, "plain", "black"),
           legend.labs = c("Elra", "Tec"),
           title = "B)",
           size = 0.5
)

PFS$plot <- PFS$plot + 
  geom_vline(xintercept = c(3, 6, 9, 12, 15), linetype = "dotted", color = "gray50") +
  annotate("text", x = 3, y = timepoint1$surv[1] - 0.08, label = paste0(signif(timepoint1$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 3, y = timepoint1$surv[2] + 0.03, label = paste0(signif(timepoint1$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 6, y = timepoint2$surv[1] - 0.08, label = paste0(signif(timepoint2$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 6, y = timepoint2$surv[2] + 0.03, label = paste0(signif(timepoint2$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 9, y = timepoint3$surv[1] - 0.08, label = paste0(signif(timepoint3$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 9, y = timepoint3$surv[2] + 0.03, label = paste0(signif(timepoint3$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 12, y = timepoint4$surv[1] - 0.08, label = paste0(signif(timepoint4$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 12, y = timepoint4$surv[2] + 0.03, label = paste0(signif(timepoint4$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 15, y = timepoint5$surv[1] - 0.08, label = paste0(signif(timepoint5$surv[1]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 15, y = timepoint5$surv[2] + 0.03, label = paste0(signif(timepoint5$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = -0.7, y = 0,
           label = paste0("Elra vs tec\n",
                          " Median: ", signif(med_time$median[1], 3), " vs ", signif(med_time$median[2], 3), " months\n",
                          " HR: ", sprintf("%.2f", HR), " (95% CI ", sprintf("%.2f", CI_lower), "-", sprintf("%.2f", CI_upper),", p=",pval, ")"),
           color = "black", hjust = 0, vjust = 0, size = 3)

PFS

km_DOR <- survfit(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ bs_ab_name, weights = weights, data = data_dor)
km_DOR
## Call: survfit(formula = Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ 
##     bs_ab_name, data = data_dor, weights = weights)
## 
##    1 observation deleted due to missingness 
##                 records    n events median 0.95LCL 0.95UCL
## bs_ab_name=elra      26 28.6   13.6   6.17     4.2      NA
## bs_ab_name=tec       51 47.5   27.4   9.90     7.0      NA
med_time <- surv_median(km_DOR)

cox_model <- coxph(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ ifelse(bs_ab_name == "elra",1,0), weights = weights, data = data_dor)
cox_summary <- summary(cox_model)

HR <- round(cox_summary$coefficients[1, "exp(coef)"], 2)
CI_lower <- round(cox_summary$conf.int[1, "lower .95"], 2)
CI_upper <- round(cox_summary$conf.int[1, "upper .95"], 2)
pval <- signif(cox_summary$coefficients[1, "Pr(>|z|)"], 2)

## Specific timepoints
timepoint1 <- summary(km_DOR,times=c(3))
timepoint2 <- summary(km_DOR,times=c(6))
timepoint3 <- summary(km_DOR,times=c(9))
timepoint4 <- summary(km_DOR,times=c(12))
timepoint5 <- summary(km_DOR,times=c(15))


DOR <- ggsurvplot(km_DOR,
           pval=FALSE,
           conf.int = FALSE,
           risk.table = TRUE, # Add risk table
           fontsize = 4,
           risk.table.col = "strata", # Change risk table color by groups
           tables.height = 0.3,
           #risk.table.fontsize = 4,
           tables.theme = theme_cleantable() +
             theme(
               axis.text.y = element_text(size = 9),
               plot.title = element_text(size = 10)
               ),
           xlab="Time (months)", ylab = "Duration of response",
           xlim = c(0,18), break.x.by = c(3),
           ylim = c(0,1), break.y.by = c(0.25),
           legend = "none",
           surv.scale="percent",
           font.main = c(10, "plain", "black"),
           font.x = c(10, "plain", "black"),
           font.y = c(10, "plain", "black"),
           font.caption = c(10, "plain", "black"),
           font.tickslab = c(10, "plain", "black"),
           legend.labs = c("Elra", "Tec"),
           title = "C)",
           size = 0.5
)

DOR$plot <- DOR$plot + 
  geom_vline(xintercept = c(3, 6, 9, 12, 15), linetype = "dotted", color = "gray50") +
  annotate("text", x = 3, y = timepoint1$surv[1] -0.08, label = paste0(signif(timepoint1$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 3, y = timepoint1$surv[2] + 0.03, label = paste0(signif(timepoint1$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 6, y = timepoint2$surv[1] - 0.08, label = paste0(signif(timepoint2$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 6, y = timepoint2$surv[2] + 0.03, label = paste0(signif(timepoint2$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 9, y = timepoint3$surv[1] - 0.08, label = paste0(signif(timepoint3$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 9, y = timepoint3$surv[2] + 0.03, label = paste0(signif(timepoint3$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 12, y = timepoint4$surv[1] - 0.08, label = paste0(signif(timepoint4$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 12, y = timepoint4$surv[2] + 0.03, label = paste0(signif(timepoint4$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 15, y = timepoint5$surv[1] - 0.08, label = paste0(signif(timepoint5$surv[1]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 15, y = timepoint5$surv[2] + 0.03, label = paste0(signif(timepoint5$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = -0.7, y = 0,
           label = paste0("Elra vs tec\n",
                          " Median: ", signif(med_time$median[1], 3), " vs ", signif(med_time$median[2], 3), " months\n",
                          " HR: ", sprintf("%.2f", HR), " (95% CI ", sprintf("%.2f", CI_lower), "-", sprintf("%.2f", CI_upper),", p=",pval, ")"),
           color = "black", hjust = 0, vjust = 0, size = 3)

DOR

6.4.4 Age ≥80 years

data1 <- IPTW.data %>%
  filter(age_at_bs_ab_initiation >=80) %>%
  transmute(
    site,
    patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate,
    bs_ab_name,
    Relapse.or.death = ifelse(progression_yes_no | death_yes_no,1,0),
    Death = ifelse(death_yes_no,1,0),
    date_of_bs_ab_step_up_d1,
    date_of_pd,
    date_of_last_contact,
    date_progression_or_death = coalesce(date_of_pd, date_of_last_contact),
    Days.to.death.or.DLC = as.numeric(date_of_last_contact - date_of_bs_ab_step_up_d1),
    Days.to.death.or.relapse.or.DLC = as.numeric(date_progression_or_death  - date_of_bs_ab_step_up_d1),
    trial_eligible = ifelse(eligible_for_corresponding_bs_ab_trial_on_first_dose_use_trial_elig_sheet_to_help_you ==1,1,0),
    weights
  )

data_dor <- IPTW.data %>%
  filter(best_response %in% c("sCR","CR","VGPR","PR")) %>%
  filter(age_at_bs_ab_initiation >=80) %>%
  transmute(
    weights,
    bs_ab_name,
    day_30_response,
    day_90_response,
    best_response,
    Relapse.or.death = ifelse(progression_yes_no | death_yes_no,1,0),
    Death = ifelse(death_yes_no,1,0),
    date_of_pd,
    date_of_last_contact,
    date_of_bs_ab_step_up_d1,
    date_progression_or_death = coalesce(date_of_pd, date_of_last_contact),
    day_30_response,
    day_90_response,
    best_response,
    time_to_best_response = case_when(
      best_response == day_30_response ~ 30,
      best_response == day_90_response ~ 90,
      .default = NA
    ),
    trial_eligible = ifelse(eligible_for_corresponding_bs_ab_trial_on_first_dose_use_trial_elig_sheet_to_help_you ==1,1,0)
  ) %>%
  filter(!is.na(time_to_best_response)) %>%
  mutate(
    Days.to.death.or.DLC = as.numeric(date_of_last_contact - date_of_bs_ab_step_up_d1) - time_to_best_response,
    Days.to.death.or.relapse.or.DLC = as.numeric(date_progression_or_death  - date_of_bs_ab_step_up_d1)- time_to_best_response,
  ) %>%
  filter(Days.to.death.or.DLC >0 & Days.to.death.or.relapse.or.DLC >0)

km_OS <- survfit(Surv(Days.to.death.or.DLC/30, Death) ~ bs_ab_name, weights = weights, data = data1)
km_OS
## Call: survfit(formula = Surv(Days.to.death.or.DLC/30, Death) ~ bs_ab_name, 
##     data = data1, weights = weights)
## 
##                 records    n events median 0.95LCL 0.95UCL
## bs_ab_name=elra      18 15.0   11.1    5.2    3.47      NA
## bs_ab_name=tec       30 32.9   10.6   17.7    9.53      NA
med_time <- surv_median(km_OS)

cox_model <- coxph(Surv(Days.to.death.or.DLC/30, Death) ~ ifelse(bs_ab_name == "elra",1,0), weights = weights, data = data1)
cox_summary <- summary(cox_model)

HR <- round(cox_summary$coefficients[1, "exp(coef)"], 2)
CI_lower <- round(cox_summary$conf.int[1, "lower .95"], 2)
CI_upper <- round(cox_summary$conf.int[1, "upper .95"], 2)
pval <- signif(cox_summary$coefficients[1, "Pr(>|z|)"], 2)

## 1-year timepoint
timepoint1 <- summary(km_OS,times=c(3))
timepoint2 <- summary(km_OS,times=c(6))
timepoint3 <- summary(km_OS,times=c(9))
timepoint4 <- summary(km_OS,times=c(12))
timepoint5 <- summary(km_OS,times=c(15))

OS <- ggsurvplot(km_OS,
           pval=FALSE,
           conf.int = FALSE,
           risk.table = TRUE, # Add risk table
           fontsize = 4,
           risk.table.col = "strata", # Change risk table color by groups
           tables.height = 0.3,
           tables.theme = theme_cleantable() +
             theme(
               axis.text.y = element_text(size = 9),
               plot.title = element_text(size = 10)
               ),
           xlab="Time (months)", ylab = "Overall survival",
           xlim = c(0,18), break.x.by = c(3),
           ylim = c(0,1), break.y.by = c(0.25),
           legend = "none",
           surv.scale="percent",
           font.main = c(10, "plain", "black"),
           font.x = c(10, "plain", "black"),
           font.y = c(10, "plain", "black"),
           font.caption = c(10, "plain", "black"),
           font.tickslab = c(10, "plain", "black"),
           legend.labs = c("Elra", "Tec"),
           title = "A)",
           size = 0.5
) 

OS$plot <- OS$plot + 
  geom_vline(xintercept = c(3, 6, 9, 12, 15), linetype = "dotted", color = "gray50") +
  annotate("text", x = 3, y = timepoint1$surv[1] - 0.1, label = paste0(signif(timepoint1$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 3, y = timepoint1$surv[2] + 0.03, label = paste0(signif(timepoint1$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 6, y = timepoint2$surv[1] - 0.08, label = paste0(signif(timepoint2$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 6, y = timepoint2$surv[2] + 0.03, label = paste0(signif(timepoint2$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 9, y = timepoint3$surv[1] - 0.08, label = paste0(signif(timepoint3$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 9, y = timepoint3$surv[2] + 0.03, label = paste0(signif(timepoint3$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 12, y = timepoint4$surv[1] - 0.08, label = paste0(signif(timepoint4$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 12, y = timepoint4$surv[2] + 0.03, label = paste0(signif(timepoint4$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 15, y = timepoint5$surv[1] - 0.08, label = paste0(signif(timepoint5$surv[1]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 15, y = timepoint5$surv[2] + 0.03, label = paste0(signif(timepoint5$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = -0.7, y = 0,
           label = paste0("Elra vs tec\n",
                          " Median: ", signif(med_time$median[1], 3), " vs ", signif(med_time$median[2], 3), " months\n",
                          " HR: ", sprintf("%.2f", HR), " (95% CI ", sprintf("%.2f", CI_lower), "-", sprintf("%.2f", CI_upper),", p=",pval, ")"),
           color = "black", hjust = 0, vjust = 0, size = 3)

OS

km_PFS <- survfit(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ bs_ab_name, weights = weights, data = data1)
km_PFS
## Call: survfit(formula = Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ 
##     bs_ab_name, data = data1, weights = weights)
## 
##                 records    n events median 0.95LCL 0.95UCL
## bs_ab_name=elra      18 15.0   11.1   3.73    1.63      NA
## bs_ab_name=tec       30 32.9   14.8   8.63    6.07      NA
med_time <- surv_median(km_PFS)

cox_model <- coxph(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ ifelse(bs_ab_name == "elra",1,0), weights = weights, data = data1)
cox_summary <- summary(cox_model)

HR <- round(cox_summary$coefficients[1, "exp(coef)"], 2)
CI_lower <- round(cox_summary$conf.int[1, "lower .95"], 2)
CI_upper <- round(cox_summary$conf.int[1, "upper .95"], 2)
pval <- signif(cox_summary$coefficients[1, "Pr(>|z|)"], 2)

## Timepoints
timepoint1 <- summary(km_PFS,times=c(3))
timepoint2 <- summary(km_PFS,times=c(6))
timepoint3 <- summary(km_PFS,times=c(9))
timepoint4 <- summary(km_PFS,times=c(12))
timepoint5 <- summary(km_PFS,times=c(15))


PFS <- ggsurvplot(km_PFS,
           pval=FALSE,
           conf.int = FALSE,
           risk.table = TRUE, # Add risk table
           fontsize = 4,
           risk.table.col = "strata", # Change risk table color by groups
           tables.height = 0.3,
           tables.theme = theme_cleantable() +
             theme(
               axis.text.y = element_text(size = 9),
               plot.title = element_text(size = 10)
               ),
           xlab="Time (months)", ylab = "Progression-free survival",
           xlim = c(0,18), break.x.by = c(3),
           ylim = c(0,1), break.y.by = c(0.25),
           legend = "none",
           surv.scale="percent",
           font.main = c(10, "plain", "black"),
           font.x = c(10, "plain", "black"),
           font.y = c(10, "plain", "black"),
           font.caption = c(10, "plain", "black"),
           font.tickslab = c(10, "plain", "black"),
           legend.labs = c("Elra", "Tec"),
           title = "B)",
           size = 0.5
)

PFS$plot <- PFS$plot + 
  geom_vline(xintercept = c(3, 6, 9, 12, 15), linetype = "dotted", color = "gray50") +
  annotate("text", x = 3, y = timepoint1$surv[1] - 0.08, label = paste0(signif(timepoint1$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 3, y = timepoint1$surv[2] + 0.03, label = paste0(signif(timepoint1$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 6, y = timepoint2$surv[1] - 0.08, label = paste0(signif(timepoint2$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 6, y = timepoint2$surv[2] + 0.03, label = paste0(signif(timepoint2$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 9, y = timepoint3$surv[1] - 0.08, label = paste0(signif(timepoint3$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 9, y = timepoint3$surv[2] + 0.03, label = paste0(signif(timepoint3$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 12, y = timepoint4$surv[1] - 0.08, label = paste0(signif(timepoint4$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 12, y = timepoint4$surv[2] + 0.03, label = paste0(signif(timepoint4$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 15, y = timepoint5$surv[1] - 0.08, label = paste0(signif(timepoint5$surv[1]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 15, y = timepoint5$surv[2] + 0.03, label = paste0(signif(timepoint5$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = -0.7, y = 0,
           label = paste0("Elra vs tec\n",
                          " Median: ", signif(med_time$median[1], 3), " vs ", signif(med_time$median[2], 3), " months\n",
                          " HR: ", sprintf("%.2f", HR), " (95% CI ", sprintf("%.2f", CI_lower), "-", sprintf("%.2f", CI_upper),", p=",pval, ")"),
           color = "black", hjust = 0, vjust = 0, size = 3)

PFS

km_DOR <- survfit(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ bs_ab_name, weights = weights, data = data_dor)
km_DOR
## Call: survfit(formula = Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ 
##     bs_ab_name, data = data_dor, weights = weights)
## 
##                 records     n events median 0.95LCL 0.95UCL
## bs_ab_name=elra       9  8.13   5.47   4.20    1.00      NA
## bs_ab_name=tec       19 20.84   7.68   6.63    5.07      NA
med_time <- surv_median(km_DOR)

cox_model <- coxph(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ ifelse(bs_ab_name == "elra",1,0), weights = weights, data = data_dor)
cox_summary <- summary(cox_model)

HR <- round(cox_summary$coefficients[1, "exp(coef)"], 2)
CI_lower <- round(cox_summary$conf.int[1, "lower .95"], 2)
CI_upper <- round(cox_summary$conf.int[1, "upper .95"], 2)
pval <- signif(cox_summary$coefficients[1, "Pr(>|z|)"], 2)

## Specific timepoints
timepoint1 <- summary(km_DOR,times=c(3))
timepoint2 <- summary(km_DOR,times=c(6))
timepoint3 <- summary(km_DOR,times=c(9))
timepoint4 <- summary(km_DOR,times=c(12))
timepoint5 <- summary(km_DOR,times=c(15))


DOR <- ggsurvplot(km_DOR,
           pval=FALSE,
           conf.int = FALSE,
           risk.table = TRUE, # Add risk table
           fontsize = 4,
           risk.table.col = "strata", # Change risk table color by groups
           tables.height = 0.3,
           #risk.table.fontsize = 4,
           tables.theme = theme_cleantable() +
             theme(
               axis.text.y = element_text(size = 9),
               plot.title = element_text(size = 10)
               ),
           xlab="Time (months)", ylab = "Duration of response",
           xlim = c(0,18), break.x.by = c(3),
           ylim = c(0,1), break.y.by = c(0.25),
           legend = "none",
           surv.scale="percent",
           font.main = c(10, "plain", "black"),
           font.x = c(10, "plain", "black"),
           font.y = c(10, "plain", "black"),
           font.caption = c(10, "plain", "black"),
           font.tickslab = c(10, "plain", "black"),
           legend.labs = c("Elra", "Tec"),
           title = "C)",
           size = 0.5
)

DOR$plot <- DOR$plot + 
  geom_vline(xintercept = c(3, 6, 9, 12, 15), linetype = "dotted", color = "gray50") +
  annotate("text", x = 3, y = timepoint1$surv[1] -0.08, label = paste0(signif(timepoint1$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 3, y = timepoint1$surv[2] + 0.03, label = paste0(signif(timepoint1$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 6, y = timepoint2$surv[1] - 0.08, label = paste0(signif(timepoint2$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 6, y = timepoint2$surv[2] + 0.03, label = paste0(signif(timepoint2$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 9, y = timepoint3$surv[1] - 0.08, label = paste0(signif(timepoint3$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 9, y = timepoint3$surv[2] + 0.03, label = paste0(signif(timepoint3$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 12, y = timepoint4$surv[1] - 0.08, label = paste0(signif(timepoint4$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 12, y = timepoint4$surv[2] + 0.03, label = paste0(signif(timepoint4$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 15, y = timepoint5$surv[1] - 0.08, label = paste0(signif(timepoint5$surv[1]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 15, y = timepoint5$surv[2] + 0.03, label = paste0(signif(timepoint5$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = -0.7, y = 0,
           label = paste0("Elra vs tec\n",
                          " Median: ", signif(med_time$median[1], 3), " vs ", signif(med_time$median[2], 3), " months\n",
                          " HR: ", sprintf("%.2f", HR), " (95% CI ", sprintf("%.2f", CI_lower), "-", sprintf("%.2f", CI_upper),", p=",pval, ")"),
           color = "black", hjust = 0, vjust = 0, size = 3)

DOR

6.5 ***FIGURE 4: PFS/OS/DOR, stratified by Age ≥80 yrs

FIGURE4 <- wrap_plots(
  (OS$plot / OS$table) + plot_layout(heights = c(4, 0.9)),
  (PFS$plot / PFS$table) + plot_layout(heights = c(4, 0.9)),
  (DOR$plot / DOR$table) + plot_layout(heights = c(4, 0.9)),
  ncol = 1  # Force vertical stacking
)

ggsave("svg/FIGURE4.svg", FIGURE4, device = "svg",
       width = 6.5,
       height = 8.5,
       units = c("in"))
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_text()`).
ggsave("jpeg/FIGURE4.jpeg", FIGURE4, device = "jpeg",
       width = 6.5,
       height = 8.5,
       dpi = 600,
       units = c("in"))
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_text()`).
ggsave("tiff/FIGURE4.tiff", FIGURE4, device = "tiff",
       width = 6.5,
       height = 8.5,
       dpi = 600,
       units = c("in"),
       compression = "lzw")
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_text()`).
knitr::include_graphics("jpeg/FIGURE4.jpeg")

7 REGRESSION

7.1 ORR, ≥CR, ≥VGPR, PFS, OS, DOR, IFS, CRS G1+, ICANS G1+, steroid use (IPTW)

options(gtsummary.add_p_footnote = FALSE, gtsummary.add_estimate_to_reference_rows = FALSE)
theme_gtsummary_compact()


preds0 <- c("Response", "BsAb", "Baseline Hgb", "Baseline LDH", "ECOG 2+", "True EMD", "Prior LOTs", "Prior BCMA")
preds1 <- c("CR/sCR", "BsAb", "Baseline Hgb", "Baseline LDH", "ECOG 2+", "True EMD", "Prior LOTs", "Prior BCMA")
preds2 <- c("BsAb", "Baseline Hgb",  "Baseline LDH", "ECOG 2+","True EMD", "Prior LOTs", "Prior BCMA")
preds3 <- c("ICANS 1+", "BsAb", "Baseline Hgb", "Baseline LDH", "ECOG 2+", "True EMD", "Prior LOTs", "Prior BCMA")
preds4 <- c("CRS 1+", "BsAb", "Baseline Hgb", "Baseline LDH", "ECOG 2+", "True EMD", "Prior LOTs", "Prior BCMA")
preds5 <- c("Steroid use", "BsAb", "Baseline Hgb", "Baseline LDH", "ECOG 2+", "True EMD", "Prior LOTs", "Prior BCMA")

preds6 <- c("VGPR", "BsAb", "Baseline Hgb", "Baseline LDH", "ECOG 2+", "True EMD", "Prior LOTs", "Prior BCMA")

data0 <- IPTW.data %>%
  transmute(
    weights,
    BsAb = factor(recode(bs_ab_name, "tec"="Tec","elra"="Elra"), levels = c("Tec","Elra")),
    HT_Risk,
    best_response,
    "Infection" = ifelse(infection == 1, 1,0),
    "CrCl <40" = ifelse(cr_cl_prior_to_step_up_dose_number_1 <40, 1,0),
    "CR/sCR" = ifelse(best_response %in% c("CR","sCR"),1,0),
    "VGPR" = ifelse(best_response %in% c("CR","sCR", "VGPR"),1,0),
    "Response" = ifelse(best_response %in% c("CR","sCR", "VGPR", "PR"),1,0),
    "True EMD" = ifelse(just_prior_to_bs_ab_initiation_is_there_true_emd_if_none_move_to_column_gk == 2, 1,0),
    "Steroid use" = ifelse(steroid_use_as_tx==1,1,0),
    "Age" = age_at_bs_ab_initiation/10,
    "ECOG" = ecog_ps_at_bs_ab_initiation_0_5,
    "ECOG 2+" = case_when(
      is.na(ecog_ps_at_bs_ab_initiation_0_5) ~ NA,
      ecog_ps_at_bs_ab_initiation_0_5 %in% c(2,3,4) ~ 1,
      TRUE ~ 0),
    "ICANS 1+" = case_when(
      is.na(max_icans_grade) ~ NA,
      max_icans_grade %in% c(1,2,3,4) ~1,
      TRUE ~ 0
      ),
    "CRS 1+" = case_when(
      is.na(max_crs_grade_0_5) ~ NA,
      max_crs_grade_0_5 %in% c(1,2,3,4) ~1,
      TRUE ~ 0
      ),
    "Prior LOTs" = number_of_prior_lines_of_therapy_before_bs_ab,
    "HRCA" = ifelse(deletion_17p_prior_to_bs_ab_0_no_1_yes | t_4_14_prior_to_bs_ab_0_no_1_yes | t_14_16_prior_to_bs_ab_0_no_1_yes | gain_amp1q21_prior_to_bs_ab_0_no_1_yes, 1, 0),
    "Prior BCMA" = prior_bcma_directed_therapy_yes_no,
    prior_bcma_directed_therapy_yes_no,
    Days.since.last.BCMA = case_when(
      is.na(prior_bcma_directed_therapy_yes_no) ~ NA,
      date_of_last_exposure_to_prior_bcma_agent > date_of_bs_ab_step_up_d1 ~ NA,
      prior_bcma_directed_therapy_yes_no == 1 ~ as.numeric(as.Date(date_of_bs_ab_step_up_d1) - as.Date(date_of_last_exposure_to_prior_bcma_agent))
    ),
    Relapse.or.death = ifelse(progression_yes_no | death_yes_no,1,0),
    Death = ifelse(death_yes_no,1,0),
    date_of_pd,
    date_of_last_contact,
    date_progression_or_death = coalesce(date_of_pd, date_of_last_contact),
    Days.to.death.or.DLC = as.numeric(date_of_last_contact - date_of_bs_ab_step_up_d1),
    Days.to.death.or.relapse.or.DLC = as.numeric(date_progression_or_death  - date_of_bs_ab_step_up_d1),
    "ALC BL" = baseline_alc,
    "CrCl BL" = cr_cl_prior_to_step_up_dose_number_1,
    "Trial ineligible" = ifelse(eligible_for_corresponding_bs_ab_trial_on_first_dose_use_trial_elig_sheet_to_help_you == 1,0,1),
    "Baseline platelets" = baseline_pl_ts,
    "Baseline Hgb" = baseline_hgb,
    "Baseline LDH" = as.double(baseline_ldh_prior_to_bs_ab_initiation)/100,
    "Baseline ferritin" = as.double(baseline_ldh_prior_to_bs_ab_initiation)/1000,
    time_to_best_response = case_when(
      best_response == day_30_response ~ 30,
      best_response == day_90_response ~ 90,
      .default = NA
      ),
    date_of_bs_ab_step_up_d1,
    bs_ab_served_only_as_pre_car_t_bridging_0_no_1_yes,
    days_to_infection_or_last_contact = as.numeric(
      pmin(as.Date(date_of_infection_74), as.Date(date_of_last_contact), na.rm = TRUE) - as.Date(date_of_bs_ab_step_up_d1)
    ),
    infection_or_death = ifelse(!is.na(date_of_infection_74) | death_yes_no == 1,1,0),
    IVIg = ifelse(ivig_use_yes_1_no_0 == 1, 1, 0),
    ALPS = ifelse(baseline_hgb >10 & baseline_ldh_prior_to_bs_ab_initiation <300,1,0),
    baseline_pl_ts,
    baseline_anc,
    baseline_hgb,
    baseline_crp_prior_to_bs_ab_initiation,
    baseline_ferritin_prior_to_bs_ab_initiation,
    CAR.HEMATOTOX = baseline_pl_ts * -0.05 + baseline_anc * -0.03 + baseline_hgb * -0.20 + as.numeric(baseline_crp_prior_to_bs_ab_initiation) * 0.15 + as.numeric(baseline_ferritin_prior_to_bs_ab_initiation) * 0.001
    ) %>%
  mutate(
    "BCMA exposure status" = factor(case_when(
      #prior_bcma_directed_therapy_yes_no == 0 ~ "  Never",
      Days.since.last.BCMA < 365 ~ "  <1 yr",
      Days.since.last.BCMA >= 365 ~ "  ≥1 yr"
    ), c("  ≥1 yr","  <1 yr")),
    `CAR-HEMATOTOX risk` = factor(
      case_when(
        CAR.HEMATOTOX < 0 ~ "Low",
        CAR.HEMATOTOX <= 5 ~ "Intermediate",
        CAR.HEMATOTOX > 5 ~ "High"
      ),
      levels = c("Low", "Intermediate", "High")
  )
  )

data1 <- data0 %>% filter(bs_ab_served_only_as_pre_car_t_bridging_0_no_1_yes==0) ## Exclude patients who received bsAb as bridging

data2 <- data1 %>% 
  filter(best_response %in% c("PR","VGPR","sCR", "CR")) %>%
  filter(!is.na(time_to_best_response)) %>%
  mutate(
    Days.to.death.or.DLC = as.numeric(date_of_last_contact - date_of_bs_ab_step_up_d1) - time_to_best_response,
    Days.to.death.or.relapse.or.DLC = as.numeric(date_progression_or_death  - date_of_bs_ab_step_up_d1)- time_to_best_response,
  ) %>%
  filter(Days.to.death.or.DLC >0 & Days.to.death.or.relapse.or.DLC >0)


## Response
uv_tab <- tbl_uvregression(
  data0 %>%
    dplyr::select(all_of(preds0)),
  method = glm,
  y = `Response`,
  method.args = list(family = binomial, weights = data0$weights),
  exponentiate = TRUE,
  add_estimate_to_reference_rows = FALSE
) %>% modify_header(estimate = "**Estimate**")

mv_tab<-glm( `Response` ~ BsAb + `True EMD` + `ECOG 2+` + `Prior LOTs` + `Prior BCMA` + `Baseline Hgb` + `Baseline LDH`, data0 %>%
    dplyr::select(all_of(preds0), weights = weights) , family = binomial ) %>%
  tbl_regression(exponentiate=TRUE) %>% modify_header(estimate = "**Estimate**") 

table_response <- tbl_merge(
  list(uv_tab, mv_tab),
  tab_spanner = c("**Univariate**", "**Multivariable**")
) 

## VGPR or better
uv_tab <- tbl_uvregression(
  data0 %>%
    dplyr::select(all_of(preds6)),
  method = glm,
  y = `VGPR`,
  method.args = list(family = binomial, weights = data0$weights),
  exponentiate = TRUE,
  add_estimate_to_reference_rows = FALSE,
) %>% modify_header(estimate = "**Estimate**")

mv_tab<-glm( `VGPR` ~ BsAb +`True EMD` + `ECOG 2+` + `Prior LOTs` + `Prior BCMA` + `Baseline Hgb` + `Baseline LDH`, data0 %>%
    dplyr::select(all_of(preds6), weights = weights) , family = binomial ) %>%
  tbl_regression(exponentiate=TRUE) %>% modify_header(estimate = "**Estimate**")

table_VGPR <- tbl_merge(
  list(uv_tab, mv_tab),
  tab_spanner = c("**Univariate**", "**Multivariable**")
) 

## sCR/CR
uv_tab <- tbl_uvregression(
  data0 %>%
    dplyr::select(all_of(preds1)),
  method = glm,
  y = `CR/sCR`,
  method.args = list(family = binomial, weights = data0$weights),
  exponentiate = TRUE,
  add_estimate_to_reference_rows = FALSE,
) %>% modify_header(estimate = "**Estimate**")

mv_tab<-glm( `CR/sCR` ~ BsAb + `True EMD` + `ECOG 2+` + `Prior LOTs` + `Prior BCMA` + `Baseline Hgb` + `Baseline LDH`, data0 %>%
    dplyr::select(all_of(preds1), weights = weights) , family = binomial ) %>%
  tbl_regression(exponentiate=TRUE) %>% modify_header(estimate = "**Estimate**")

table_sCR.CR <- tbl_merge(
  list(uv_tab, mv_tab),
  tab_spanner = c("**Univariate**", "**Multivariable**")
) 

## OS

uv_tab <- tbl_uvregression(
  data1[c(preds2)],
  method = coxph,
  y = Surv(data1$Days.to.death.or.DLC, data1$Death),
  exponentiate = TRUE,
  method.args = list(weights = data1$weights)
) %>% modify_header(estimate = "**Estimate**")

mv_tab<-coxph( Surv(Days.to.death.or.DLC, Death) ~ BsAb + `True EMD` + `ECOG 2+` + `Prior LOTs` + `Prior BCMA` + `Baseline Hgb` + `Baseline LDH`, data1,  weights = weights) %>%
  tbl_regression(exponentiate=TRUE) %>% modify_header(estimate = "**Estimate**")

table_os <- tbl_merge(
  list(uv_tab, mv_tab),
  tab_spanner = c("**Univariate**", "**Multivariable**")
) 


## PFS
uv_tab <- tbl_uvregression(
  data1[c(preds2)],
  method = coxph,
  y = Surv(data1$Days.to.death.or.relapse.or.DLC, data1$Relapse.or.death),
  exponentiate = TRUE,
  method.args = list(weights = data1$weights)
) %>% modify_header(estimate = "**Estimate**")

mv_tab<-coxph( Surv(Days.to.death.or.relapse.or.DLC, Relapse.or.death) ~ BsAb + `True EMD` + `ECOG 2+` + `Prior LOTs` + `Prior BCMA` + `Baseline Hgb` + `Baseline LDH`, weights = weights, data1) %>%
  tbl_regression(exponentiate=TRUE) %>% modify_header(estimate = "**Estimate**")

table_pfs <- tbl_merge(
  list(uv_tab, mv_tab),
  tab_spanner = c("**Univariate**", "**Multivariable**")
)

## DOR
uv_tab <- tbl_uvregression(
  data2[c(preds2)],
  method = coxph,
  y = Surv(data2$Days.to.death.or.relapse.or.DLC, data2$Relapse.or.death),
  exponentiate = TRUE,
  method.args = list(weights = data2$weights)
) %>% modify_header(estimate = "**Estimate**")

mv_tab<-coxph( Surv(  Days.to.death.or.relapse.or.DLC, Relapse.or.death) ~ BsAb + `True EMD` + `ECOG 2+` + `Prior LOTs` + `Prior BCMA` + `Baseline Hgb` + `Baseline LDH`, weights = weights, data2) %>%
  tbl_regression(exponentiate=TRUE) %>% modify_header(estimate = "**Estimate**")

table_dor <- tbl_merge(
  list(uv_tab, mv_tab),
  tab_spanner = c("**Univariate**", "**Multivariable**")
)

## IFS
uv_tab <- tbl_uvregression(
  data1[c(preds2)],
  method = coxph,
  y = Surv(data1$days_to_infection_or_last_contact, data1$infection_or_death),
  exponentiate = TRUE,
  method.args = list(weights = data1$weights)
 ) %>% modify_header(estimate = "**Estimate**")

mv_tab<-coxph( Surv(  days_to_infection_or_last_contact, infection_or_death) ~ BsAb + `True EMD` + `ECOG 2+` + `Prior LOTs` + `Prior BCMA` + `Baseline Hgb` + `Baseline LDH`, weights = weights, data1) %>%
  tbl_regression(exponentiate=TRUE) %>% modify_header(estimate = "**Estimate**")

table_ifs <- tbl_merge(
  list(uv_tab, mv_tab),
  tab_spanner = c("**Univariate**", "**Multivariable**")
)

## ICANS 1+
uv_tab <- tbl_uvregression(
  data0 %>%
    dplyr::select(all_of(preds3)),
  method = glm,
  y = `ICANS 1+`,
  method.args = list(family = binomial, weights = data0$weights),
  exponentiate = TRUE,
  add_estimate_to_reference_rows = FALSE
) %>% modify_header(estimate = "**Estimate**")

mv_tab<-glm( `ICANS 1+` ~ BsAb + `True EMD` + `ECOG 2+` + `Prior LOTs` + `Prior BCMA` + `Baseline Hgb` + `Baseline LDH`, data0 %>%
    dplyr::select(all_of(preds3), weights = weights) , family = binomial ) %>%
  tbl_regression(exponentiate=TRUE) %>% modify_header(estimate = "**Estimate**")

table_ICANS <- tbl_merge(
  list(uv_tab, mv_tab),
  tab_spanner = c("**Univariate**", "**Multivariable**")
) 

## CRS 1+
uv_tab <- tbl_uvregression(
  data0 %>%
    dplyr::select(all_of(preds4)),
  method = glm,
  y = `CRS 1+`,
  method.args = list(family = binomial, weights = data0$weights),
  exponentiate = TRUE,
  add_estimate_to_reference_rows = FALSE
) %>% modify_header(estimate = "**Estimate**")

mv_tab<-glm( `CRS 1+` ~ BsAb + `True EMD` + `ECOG 2+` + `Prior LOTs` + `Prior BCMA` + `Baseline Hgb` + `Baseline LDH`, data0 %>%
    dplyr::select(all_of(preds4), weights = weights) , family = binomial ) %>%
  tbl_regression(exponentiate=TRUE) %>% modify_header(estimate = "**Estimate**")

table_CRS <- tbl_merge(
  list(uv_tab, mv_tab),
  tab_spanner = c("**Univariate**", "**Multivariable**")
) 

## Steroid use
uv_tab <- tbl_uvregression(
  data0 %>%
    dplyr::select(all_of(preds5)),
  method = glm,
  y = `Steroid use`,
  method.args = list(family = binomial, weights = data0$weights),
  exponentiate = TRUE,
  add_estimate_to_reference_rows = FALSE
) %>% modify_header(estimate = "**Estimate**")

mv_tab<-glm( `Steroid use` ~ BsAb + `True EMD` + `ECOG 2+` + `Prior LOTs` + `Prior BCMA` + `Baseline Hgb` + `Baseline LDH`, data0 %>%
    dplyr::select(all_of(preds5), weights = weights) , family = binomial ) %>%
  tbl_regression(exponentiate=TRUE) %>% modify_header(estimate = "**Estimate**")

table_steroid <- tbl_merge(
  list(uv_tab, mv_tab),
  tab_spanner = c("**Univariate**", "**Multivariable**")
) 


# Combine all tables into one
final_table <- tbl_stack(
  list(
    table_response %>% modify_header() %>% modify_footnote(everything() ~ NA, abbreviation = TRUE) %>%
      modify_footnote(estimate_1 ~ "Estimate = Odds ratio or hazard ratio", abbreviation = TRUE),
    table_VGPR,
    table_sCR.CR,
    table_os,
    table_pfs,
    table_dor,
    table_ifs,
    table_ICANS,
    table_CRS,
    table_steroid
  ),
  group_header = c(
    "Response",
    "≥VGPR",
    "CR/sCR",
    "Overall Survival",
    "Progression-free Survival",
    "Duration of Response",
    "Infection-free Survival",
    "ICANS 1+",
    "CRS 1+",
    "Steroid use"
  )
)
  
final_table
Characteristic
Univariate
Multivariable
N Estimate1 95% CI p-value Estimate 95% CI p-value
Response
BsAb 370





    Tec
— —
— —
    Elra
0.86 0.55, 1.35 0.5 0.86 0.52, 1.42 0.6
Baseline Hgb 370 1.46 1.29, 1.65 <0.001 1.34 1.18, 1.53 <0.001
Baseline LDH 370 0.86 0.77, 0.94 0.002 0.86 0.77, 0.95 0.004
ECOG 2+ 370 0.39 0.25, 0.61 <0.001 0.42 0.26, 0.69 <0.001
True EMD 370 0.91 0.52, 1.63 0.8 0.94 0.50, 1.78 0.9
Prior LOTs 370 0.87 0.79, 0.94 0.001 0.91 0.82, 1.01 0.087
Prior BCMA 370 0.56 0.36, 0.85 0.008 0.50 0.29, 0.87 0.014
≥VGPR
BsAb 370





    Tec
— —
— —
    Elra
0.62 0.40, 0.95 0.031 0.57 0.35, 0.93 0.026
Baseline Hgb 370 1.38 1.23, 1.55 <0.001 1.32 1.17, 1.49 <0.001
Baseline LDH 370 0.84 0.73, 0.93 0.003 0.84 0.73, 0.93 0.004
ECOG 2+ 370 0.52 0.33, 0.81 0.004 0.57 0.34, 0.93 0.026
True EMD 370 0.63 0.35, 1.10 0.11 0.56 0.30, 1.04 0.071
Prior LOTs 370 0.88 0.80, 0.95 0.003 0.95 0.86, 1.06 0.4
Prior BCMA 370 0.47 0.31, 0.72 <0.001 0.41 0.24, 0.70 0.001
CR/sCR
BsAb 370





    Tec
— —
— —
    Elra
0.70 0.44, 1.10 0.12 0.68 0.40, 1.12 0.13
Baseline Hgb 370 1.44 1.29, 1.63 <0.001 1.39 1.23, 1.58 <0.001
Baseline LDH 370 0.86 0.74, 0.96 0.019 0.87 0.75, 0.98 0.044
ECOG 2+ 370 0.51 0.31, 0.81 0.006 0.61 0.36, 1.03 0.068
True EMD 370 0.66 0.35, 1.20 0.2 0.62 0.31, 1.19 0.2
Prior LOTs 370 0.89 0.81, 0.97 0.014 0.95 0.85, 1.06 0.4
Prior BCMA 370 0.56 0.37, 0.86 0.009 0.52 0.31, 0.89 0.018
Overall Survival
BsAb 352





    Tec
— —
— —
    Elra
1.41 0.96, 2.07 0.083 1.46 0.97, 2.20 0.073
Baseline Hgb 352 0.73 0.66, 0.81 <0.001 0.77 0.69, 0.86 <0.001
Baseline LDH 352 1.13 1.06, 1.20 <0.001 1.15 1.10, 1.20 <0.001
ECOG 2+ 352 2.25 1.61, 3.16 <0.001 1.92 1.32, 2.78 <0.001
True EMD 352 1.34 0.88, 2.03 0.2 1.28 0.82, 2.00 0.3
Prior LOTs 352 1.06 0.99, 1.13 0.091 1.04 0.97, 1.11 0.3
Prior BCMA 352 1.09 0.77, 1.54 0.6 1.20 0.81, 1.78 0.4
Progression-free Survival
BsAb 355





    Tec
— —
— —
    Elra
1.43 1.05, 1.94 0.021 1.52 1.11, 2.09 0.009
Baseline Hgb 355 0.77 0.71, 0.83 <0.001 0.81 0.74, 0.89 <0.001
Baseline LDH 355 1.18 1.08, 1.28 <0.001 1.17 1.08, 1.27 <0.001
ECOG 2+ 355 1.63 1.23, 2.16 <0.001 1.43 1.04, 1.95 0.027
True EMD 355 1.04 0.74, 1.47 0.8 0.94 0.67, 1.32 0.7
Prior LOTs 355 1.10 1.05, 1.16 <0.001 1.06 1.00, 1.12 0.056
Prior BCMA 355 1.58 1.18, 2.10 0.002 1.63 1.17, 2.27 0.004
Duration of Response
BsAb 177





    Tec
— —
— —
    Elra
1.49 0.87, 2.56 0.14 1.50 0.85, 2.63 0.2
Baseline Hgb 177 0.84 0.74, 0.95 0.007 0.83 0.72, 0.96 0.014
Baseline LDH 177 1.05 0.90, 1.23 0.5 1.04 0.89, 1.21 0.7
ECOG 2+ 177 1.05 0.64, 1.72 0.8 0.89 0.52, 1.52 0.7
True EMD 177 1.48 0.75, 2.93 0.3 1.36 0.68, 2.72 0.4
Prior LOTs 177 1.01 0.89, 1.15 0.8 1.02 0.90, 1.17 0.7
Prior BCMA 177 1.08 0.67, 1.74 0.8 1.06 0.65, 1.75 0.8
Infection-free Survival
BsAb 358





    Tec
— —
— —
    Elra
1.25 0.93, 1.67 0.13 1.24 0.92, 1.67 0.2
Baseline Hgb 358 0.87 0.81, 0.94 <0.001 0.90 0.83, 0.97 0.010
Baseline LDH 358 1.10 1.05, 1.15 <0.001 1.09 1.05, 1.14 <0.001
ECOG 2+ 358 1.70 1.30, 2.22 <0.001 1.52 1.14, 2.02 0.004
True EMD 358 1.16 0.83, 1.62 0.4 1.10 0.79, 1.54 0.6
Prior LOTs 358 0.98 0.93, 1.03 0.4 0.98 0.93, 1.04 0.5
Prior BCMA 358 0.85 0.65, 1.11 0.2 0.99 0.73, 1.35 >0.9
ICANS 1+
BsAb 370





    Tec
— —
— —
    Elra
0.87 0.49, 1.52 0.6 0.91 0.49, 1.65 0.8
Baseline Hgb 370 0.83 0.72, 0.95 0.008 0.84 0.72, 0.98 0.024
Baseline LDH 370 1.05 0.95, 1.15 0.3 1.04 0.93, 1.14 0.5
ECOG 2+ 370 3.03 1.77, 5.21 <0.001 2.44 1.39, 4.31 0.002
True EMD 370 1.40 0.69, 2.69 0.3 1.31 0.62, 2.62 0.5
Prior LOTs 370 0.88 0.77, 0.98 0.034 0.88 0.76, 1.01 0.069
Prior BCMA 370 0.61 0.36, 1.03 0.067 1.01 0.53, 1.92 >0.9
CRS 1+
BsAb 370





    Tec
— —
— —
    Elra
0.53 0.34, 0.82 0.004 0.55 0.35, 0.86 0.010
Baseline Hgb 370 0.94 0.85, 1.04 0.3 0.91 0.81, 1.01 0.076
Baseline LDH 370 0.98 0.90, 1.06 0.5 0.97 0.89, 1.06 0.5
ECOG 2+ 370 1.08 0.70, 1.68 0.7 0.97 0.60, 1.54 0.9
True EMD 370 0.58 0.33, 1.02 0.061 0.53 0.29, 0.94 0.033
Prior LOTs 370 0.92 0.84, 1.00 0.045 0.87 0.79, 0.96 0.006
Prior BCMA 370 1.17 0.77, 1.76 0.5 1.84 1.12, 3.04 0.016
Steroid use
BsAb 355





    Tec
— —
— —
    Elra
0.83 0.48, 1.40 0.5 0.87 0.50, 1.50 0.6
Baseline Hgb 355 0.89 0.79, 1.01 0.074 0.90 0.79, 1.03 0.14
Baseline LDH 355 1.07 0.98, 1.17 0.12 1.05 0.96, 1.15 0.3
ECOG 2+ 355 1.71 1.03, 2.83 0.036 1.32 0.78, 2.23 0.3
True EMD 355 1.61 0.85, 2.97 0.13 1.57 0.81, 2.94 0.2
Prior LOTs 355 0.88 0.79, 0.98 0.028 0.94 0.83, 1.06 0.3
Prior BCMA 355 0.47 0.29, 0.78 0.003 0.65 0.36, 1.16 0.15
1 Estimate = Odds ratio or hazard ratio

8 NON-LINEAR ANALYSIS

8.1 OS based on age

data1 <- IPTW.data %>%
  transmute(
    site,
    age_at_bs_ab_initiation,
    patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate,
    bs_ab_name = factor(recode(bs_ab_name, "tec"="Tec","elra"="Elra"), levels = c("Tec","Elra")),
    Relapse.or.death = ifelse(progression_yes_no | death_yes_no,1,0),
    Death = ifelse(death_yes_no,1,0),
    date_of_bs_ab_step_up_d1,
    date_of_pd,
    date_of_last_contact,
    date_progression_or_death = coalesce(date_of_pd, date_of_last_contact),
    Days.to.death.or.DLC = as.numeric(date_of_last_contact - date_of_bs_ab_step_up_d1),
    Days.to.death.or.relapse.or.DLC = as.numeric(date_progression_or_death  - date_of_bs_ab_step_up_d1),
    trial_eligible = ifelse(eligible_for_corresponding_bs_ab_trial_on_first_dose_use_trial_elig_sheet_to_help_you ==1,1,0),
    weights
  )


dd <- datadist(data1)
options(datadist='dd')

ref_value = 85

S <- Surv(data1$Days.to.death.or.DLC, data1$Death)
f <- cph(S ~ rcs(age_at_bs_ab_initiation,3) * bs_ab_name , x=TRUE, y=TRUE,surv=TRUE,data=data1, weights = weights)

coxph( S ~ age_at_bs_ab_initiation * bs_ab_name , data=data1) %>% tbl_regression(exponentiate=TRUE)
Characteristic HR1 95% CI1 p-value
age_at_bs_ab_initiation 0.99 0.97, 1.01 0.3
bs_ab_name


    Tec — —
    Elra 0.05 0.00, 1.01 0.051
age_at_bs_ab_initiation * bs_ab_name


    age_at_bs_ab_initiation * Elra 1.05 1.00, 1.09 0.029
1 HR = Hazard Ratio, CI = Confidence Interval
model <- Predict(f, bs_ab_name, age_at_bs_ab_initiation, fun=function(x) exp(x) /exp(predict(f, data.frame(age_at_bs_ab_initiation = ref_value, bs_ab_name = "Tec")) )    )

OS_age <- ggplot(as.data.frame(model),aes(x=age_at_bs_ab_initiation, y=yhat, z = bs_ab_name)) + 
  geom_line(
    aes(col = bs_ab_name)
  ) + 
  geom_ribbon(data = filter(model, bs_ab_name =="Elra"), aes(ymin=lower, ymax=upper), alpha=0.2, linetype=0, fill = "steelblue" ) + 
  geom_ribbon(data = filter(model, bs_ab_name =="Tec"), aes(ymin=lower, ymax=upper), alpha=0.2, linetype=0, fill = "coral2" ) + 
  theme_bw() + 
  labs(
    title = "A)",
    x = "Age (years)",
    y = "OS (HR)") +
  scale_color_manual(values = c ("coral2", "steelblue")) + 
  theme(title = element_text(size=16,face="bold"),legend.title=element_blank()) + 
  #coord_cartesian(ylim = c(0.5, 5.5)) +
  geom_hline(yintercept = 1, linetype = "dashed", color = "gray", size = 0.5)
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
OS_age

8.2 PFS based on age

data1 <- IPTW.data %>%
  transmute(
    site,
    age_at_bs_ab_initiation,
    patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate,
    bs_ab_name = factor(recode(bs_ab_name, "tec"="Tec","elra"="Elra"), levels = c("Tec","Elra")),
    Relapse.or.death = ifelse(progression_yes_no | death_yes_no,1,0),
    Death = ifelse(death_yes_no,1,0),
    date_of_bs_ab_step_up_d1,
    date_of_pd,
    date_of_last_contact,
    date_progression_or_death = coalesce(date_of_pd, date_of_last_contact),
    Days.to.death.or.DLC = as.numeric(date_of_last_contact - date_of_bs_ab_step_up_d1),
    Days.to.death.or.relapse.or.DLC = as.numeric(date_progression_or_death  - date_of_bs_ab_step_up_d1),
    trial_eligible = ifelse(eligible_for_corresponding_bs_ab_trial_on_first_dose_use_trial_elig_sheet_to_help_you ==1,1,0),
    weights
  )


dd <- datadist(data1)
options(datadist='dd')

ref_value = 85

S <- Surv(data1$Days.to.death.or.relapse.or.DLC, data1$Relapse.or.death)
f <- cph(S ~ rcs(age_at_bs_ab_initiation,3) * bs_ab_name , x=TRUE, y=TRUE,surv=TRUE,data=data1, weights = weights)

coxph( S ~ age_at_bs_ab_initiation * bs_ab_name , data=data1) %>% tbl_regression(exponentiate=TRUE)
Characteristic HR1 95% CI1 p-value
age_at_bs_ab_initiation 0.98 0.97, 1.00 0.034
bs_ab_name


    Tec — —
    Elra 0.54 0.06, 4.75 0.6
age_at_bs_ab_initiation * bs_ab_name


    age_at_bs_ab_initiation * Elra 1.01 0.98, 1.04 0.4
1 HR = Hazard Ratio, CI = Confidence Interval
model <- Predict(f, bs_ab_name, age_at_bs_ab_initiation, fun=function(x) exp(x) /exp(predict(f, data.frame(age_at_bs_ab_initiation = ref_value, bs_ab_name = "Tec")) )    )

PFS_age <- ggplot(as.data.frame(model),aes(x=age_at_bs_ab_initiation, y=yhat, z = bs_ab_name)) + 
  geom_line(
    aes(col = bs_ab_name)
  ) + 
  geom_ribbon(data = filter(model, bs_ab_name =="Elra"), aes(ymin=lower, ymax=upper), alpha=0.2, linetype=0, fill = "steelblue" ) + 
  geom_ribbon(data = filter(model, bs_ab_name =="Tec"), aes(ymin=lower, ymax=upper), alpha=0.2, linetype=0, fill = "coral2" ) + 
  theme_bw() + 
  labs(
    title = "B)",
    x = "Age (years)",
    y = "PFS (HR)") +
  scale_color_manual(values = c ("coral2", "steelblue")) + 
  theme(title = element_text(size=16,face="bold"),legend.title=element_blank()) + 
  #coord_cartesian(ylim = c(0.5, 5.5)) +
  geom_hline(yintercept = 1, linetype = "dashed", color = "gray", size = 0.5)

PFS_age

8.3 ***FIGURE 5: OS/PFS and age

combined_age <- (OS_age / PFS_age) +
  plot_layout(guides = "collect", axes = "collect")


ggsave("svg/FIGURE5.svg", combined_age, device = "svg",
       width = 6.5,
       height = 6.5,
       units = c("in"))

ggsave("jpeg/FIGURE5.jpeg", combined_age, device = "jpeg",
       width = 6.5,
       height = 6.5,
       dpi = 600,
       units = c("in"))

ggsave("tiff/FIGURE5.tiff", combined_age, device = "tiff",
       width = 6.5,
       height = 6.5,
       dpi = 600,
       units = c("in"),
       compression = "lzw")

knitr::include_graphics("jpeg/FIGURE5.jpeg")

9 CONTOUR PLOTS

9.1 FIGURE 6: Baseline Hgb and LDH

data1 <- IPTW.data %>%
  filter(bs_ab_served_only_as_pre_car_t_bridging_0_no_1_yes==0) %>% ## Exclude patients who received bsAb as bridging
  #filter(bs_ab_name == "tec") %>%
  transmute(
    bs_ab_name,
    weights,
    "CR/sCR" = ifelse(best_response %in% c("CR", "sCR"),1,0),
    "PD" = ifelse(best_response %in% c("PD"),1,0),
    "True EMD" = ifelse(just_prior_to_bs_ab_initiation_is_there_true_emd_if_none_move_to_column_gk == 2, 1,0),
    "Age" = age_at_bs_ab_initiation/10,
    "ECOG" = ecog_ps_at_bs_ab_initiation_0_5,
    ECOG_2_4 = case_when(
      is.na(ecog_ps_at_bs_ab_initiation_0_5) ~ NA,
      ecog_ps_at_bs_ab_initiation_0_5 %in% c(2,3,4) ~ 1,
      TRUE ~ 0),
    "Prior LOTs" = number_of_prior_lines_of_therapy_before_bs_ab,
    "HRCA" = ifelse(deletion_17p_prior_to_bs_ab_0_no_1_yes | t_4_14_prior_to_bs_ab_0_no_1_yes | t_14_16_prior_to_bs_ab_0_no_1_yes | gain_amp1q21_prior_to_bs_ab_0_no_1_yes, 1, 0),
    "Prior BCMA" = prior_bcma_directed_therapy_yes_no,
    Relapse.or.death = ifelse(progression_yes_no | death_yes_no,1,0),
    Death = ifelse(death_yes_no,1,0),
    date_of_pd,
    date_of_last_contact,
    date_progression_or_death = coalesce(date_of_pd, date_of_last_contact),
    Days.to.death.or.DLC = as.numeric(date_of_last_contact - date_of_bs_ab_step_up_d1),
    Days.to.death.or.relapse.or.DLC = as.numeric(date_progression_or_death  - date_of_bs_ab_step_up_d1),
    "ALC BL" = baseline_alc,
    "CrCl BL" = cr_cl_prior_to_step_up_dose_number_1,
    "Trial ineligible" = ifelse(eligible_for_corresponding_bs_ab_trial_on_first_dose_use_trial_elig_sheet_to_help_you == 1,0,1),
    "Baseline platelets" = baseline_pl_ts,
    baseline_hgb,
    baseline_ldh = as.double(baseline_ldh_prior_to_bs_ab_initiation),
    baseline_ferritin = as.numeric(baseline_ferritin_prior_to_bs_ab_initiation),
    ) %>% filter(!is.na(baseline_ldh))


## Tec
data_tec = data1 %>% filter(bs_ab_name == "tec")

dd_tec <- datadist(data_tec)
options(datadist='dd_tec')

S_tec <- Surv(data_tec$Days.to.death.or.relapse.or.DLC/30, data_tec$Relapse.or.death)
f_tec <- cph(S_tec ~ rcs(baseline_hgb, 3) * rcs(baseline_ldh, 3) , x=TRUE, y=TRUE,surv=TRUE,data=data_tec, weights=weights)
med_tec <- Quantile(f_tec)

pred_tec  <- data.frame(Predict(f_tec,  baseline_hgb, baseline_ldh,
                                fun = function(x) med_tec(lp = x))) %>%
  mutate(drug = "Teclistamab")

## Elra
data_elra = data1 %>% filter(bs_ab_name == "elra")

dd_elra <- datadist(data_elra)
options(datadist='dd_elra')

S_elra <- Surv(data_elra$Days.to.death.or.relapse.or.DLC/30, data_elra$Relapse.or.death)
f_elra <- cph(S_elra ~ rcs(baseline_hgb, 3) * rcs(baseline_ldh, 3) , x=TRUE, y=TRUE,surv=TRUE,data=data_elra, weights=weights)
med_elra <- Quantile(f_elra)

pred_elra <- data.frame(Predict(f_elra, baseline_hgb, baseline_ldh,
                                fun = function(x) med_elra(lp = x))) %>%
  mutate(drug = "Elranatamab")

## Combine both predictions
all_pred <- bind_rows(pred_tec, pred_elra)
z_range <- range(all_pred$yhat, na.rm = TRUE)
brks    <- pretty(z_range, n = 10)   # tweak n if you want more/fewer bands

all_pred <- all_pred %>%
  mutate(
    z_band = cut(
      yhat,
      breaks        = brks,
      include.lowest = TRUE,
      right          = FALSE
    )
  )

z_levels <- levels(all_pred$z_band)

cols <- viridis_pal()(length(z_levels))


contour_plot_tec <- all_pred %>%
  filter(drug == "Teclistamab") %>%
  ggplot(aes(baseline_hgb, baseline_ldh, fill = z_band)) +
  geom_raster() +   # or geom_tile(); blocky but guarantees alignment
  theme_classic() +
  theme(
    panel.grid.major = element_line(color = "grey85"),
    panel.grid.minor = element_line(color = "grey92")
  )+
  labs(
    title = "A) Teclistamab",
    x = "Baseline Hgb (g/dL)",
    y = "Baseline LDH (U/L)",
    fill = "Median PFS (months)"
  ) +
  scale_fill_manual(
    values       = cols,
    breaks       = z_levels,
    drop         = FALSE,
    na.translate = FALSE
  ) +
  scale_x_continuous(breaks = seq(0, 20, 1))+
  scale_y_continuous(breaks = seq(0, 2000, 100))+
  coord_cartesian(xlim = c(7, 13), ylim = c(120, 780))

contour_plot_tec

contour_plot_elra <- all_pred %>%
  filter(drug == "Elranatamab") %>%
  ggplot(aes(baseline_hgb, baseline_ldh, fill = z_band)) +
  geom_raster() +
  theme_classic() +
  theme(
    legend.position = "none",
    panel.grid.major = element_line(color = "grey85"),
    panel.grid.minor = element_line(color = "grey92")
    ) +
  labs(
    title = "B) Elranatamab",
    x = "Baseline Hgb (g/dL)",
    y = "Baseline LDH (U/L)"
    #fill = "Median PFS (months)"
  ) +
  scale_fill_manual(
    values       = cols,
    breaks       = z_levels,
    drop         = FALSE,
    na.translate = FALSE
  ) +
  scale_x_continuous(breaks = seq(0, 20, 1))+
  scale_y_continuous(breaks = seq(0, 2000, 100))+
  coord_cartesian(xlim = c(7, 13), ylim = c(120, 780))

 contour_plot_elra

9.2 ***FIGURE 6: Hgb and LDH contour

combined_contour <- (contour_plot_tec / contour_plot_elra) +
  plot_layout(guides = "collect", axes = "collect")


ggsave("svg/FIGURE6.svg", combined_contour, device = "svg",
       width = 6.5,
       height = 6.5,
       units = c("in"))

ggsave("jpeg/FIGURE6.jpeg", combined_contour, device = "jpeg",
       width = 6.5,
       height = 6.5,
       dpi = 600,
       units = c("in"))

ggsave("tiff/FIGURE6.tiff", combined_contour, device = "tiff",
       width = 6.5,
       height = 6.5,
       dpi = 600,
       units = c("in"),
       compression = "lzw")

knitr::include_graphics("jpeg/FIGURE6.jpeg")